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ABSTRACT 

A  computer  program  to  perform  dynamic  systems  analyses  and 
plot  the  results  is  presented.  Linear  systems  and  a  large  classifica¬ 
tion  of  nonlinear  systems  representing  engineering,  scientific,  and 
economic  disciplines  can  be  modeled  to  permit  application  of  the  computer 
program.  Two  examples  are  given  to  demonstrate  the  capabilities  of  the 
analysis  tool.  The  mathematical  model  of  a  missile  guidance  and  control 
system  is  analyzed  and  a  ratio  of  polynomials  representing  the  closed 
loop  transfer  function  of  a  high  performance  model  follower  aircraft  is 
evaluated.  Linear  differential  equations  to  the  100th  order  having 
real  or  complex  roots  can  be  studied.  System  characteristic  equation 
root  loci  and  system  transfer  functions  are  plotted. 
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INTRODUCTION 


Systems  analyses  require  that  mathematical  models  be  developed 
and  exercised  to  determine  cause-and-effect  relationships.  One  powerful 
method  is  to  generate  linearized  equations  for  each  physical  phenomenon 
occurring  in  the  process  and  combine  these  into  input-output  relation¬ 
ships  through  the  application  of  Laplace  transformations. 

Studying  the  characteristic  equations  of  the  input-output 
relations  determines  parametric  stability  trends. Steady  state  fre¬ 
quency  response  of  the  system  provides  another  useful  characterization 
of  system  stability.  Application  of  the  inverse  Laplace  transformation 

to  the  process  transfer  functions  provides  transient  histories  of  output 

(2) 

for  specific  inputs. 

This  report  discusses  a  program  that  mechanizes  another  part 
of  the  analytical  procedure,  the  determination  of  the  frequency  response 
of  linear  systems  and  the  modified  phase-amplitude  characteristic  for 
nonlinear  systems.  Reference  1,  published  in  November  1966,  discusses 
a  program  that  mechanized  the  determination  of  the  roots  of  the  char¬ 
acteristic  equation.  Reference  2,  discusses  a  program  that  determines 
the  transient  response  of  the  process  when  it  is  subjected  to  input 
forcing  functions. 

Application  of  this  computer  program,  which  is  an  extension 
of  the  material  contained  in  Reference  1,  along  with  the  computer  pro¬ 
gram  presented  in  Reference  2,  permits  one  to  determine  quickly  and 
accurately  the  stability  and  performance  characteristics  of  high-order 
linear  dynamic  processes  and  to  qualitatively  predict  the  stability 
and  performance  of  a  large  class  of  nonlinear  systems. 

The  techniques  developed  in  this  report  have  been  applied  to 
specific  studies. A  more  extensive  utilization  of  the  methods 
will  be  contained  in  a  forthcoming  report  which  expands  Chapter  4.2  of 
Reference  5»  Preliminary  work  has  been  commenced  to  explore  the  utility 
of  these  methods  in  providing  solutions  to  optimal  inventory 
problems. 
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A  useful  linear  system  analysis  method  is  to  determine  the 
steady  state  frequency  response  of  the  system.  Evaluation  of  the  poly¬ 
nomials  representing  the  system’s  transfer  function  is  a  time-consuming 
task.  Conversion  of  the  tabulated  results  to  a  graphical  display  is 
then  required  before  conclusions  relevant  to  system  stability  can  be 
appreciated. 

The  analysis  of  nonlinear  systems  is  a  difficult  task,  but 
for  certain  classes  of  systems  it  is  possible  to  obtain  qualitative 
indications  of  stability.  Similarly,  the  task  of  obtaining  the  tabular 
data  and  converting  it  to  a  graphical  form  before  conclusions  can  be 
drawn  is  a  laborious  process. 

These  two  methods  (one  linear  and  one  nonlinear)  have  been 
coded  in  a  FORTRAN  IV  program  and  integrated  into  the  root  locus  methods 
discussed  in  Reference  1.  Main  features  of  the  program  are: 

1.  FORTRAN  IV  program.  No  machine-oriented  or  object  language. 

2.  No  FORTRAN  complex  type  statements  necessary. 

3.  No  special  graphical  plotting  equipment. 

k.  Order  of  polynomial  may  be  up  to  100. 

5-  Number  of  variations  of  given  parameter  may  be  up  to  100. 

Main  features  of  the  graphical  displays  are: 

ROOT  LOCUS  Method^ 

l.  Log  plot  of  the  third  and  fourth  quadrants  of  the  complex 
frequency  plane  from  0  to  10,000  radians/second. 

2.  Linear  plots  of  selected  regions  of  the  third  and  fourth 
quadrants  of  the  complex  frequency  plane  with  arbitrary  scales. 

LINEAR  FREQUENCY  RESPONSE  Method 

1.  Plot  of  Magnitude  (in  decibels)  and  phase  (in  degrees)  of 
transfer  function  vs  6  decades  of  frequency,  beginning  at  a  selectable 
minimum  frequency. 

2.  Linear  frequency  plots  of  magnitude  (in  decibels)  and 
phase  (in  degrees)  of  transfer  function  for  selectable  ranges  of  frequency. 
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NONLINEAR  FREQUENCY  RESPONSE  Method 

1.  Plot  of  ’Modified  Magnitude  Characteristic"  decibels) 
and  "Modified  Phase  Characteristic"  (in  degrees)  of  transfer  function  vs 
€  decades  of  frequency,  beginning  at  same  selectable  minimum  frequency 
as  linear  frequency  response  method. 

2.  Linear  frequency  plots  of  ’Modified  Magnitude  Charac¬ 
teristics"  (in  decibels)  and  "Modified  Phase  Characteristic"  (in  degrees) 
of  transfer  function  for  same  selectable  ranges  of  frequency  as  linear 
frequency  response  methods. 


2.  THE  PROBLEM 


2.1  Linear  System. 

Regardless  of  the  complexity  of  a  linear  closed- loop  system, 
its  transfer  function  can  be  reduced  to  the  equivalent  fonn  shown  in 
Figure  1.  For  multiple  loop  systems,  the  G*s  and  H’s  are  readily  ex¬ 
pressed  as  sums  of  products  of  polynomials  which  are  identified  with 
individual  elements  making  up  the  complete  system. 


KpG|  (S) 

G2(S) 

r 

F 

Kc  H,  (S) 

H2(S) 

Figure  1. 

Linear  Closed  Loop  System 

R  = 

system  input, 

C  = 

system  output, 

F  = 

system  feedback, 

E  = 

system  error. 

K  = 
p 

process  gain,  and 

K  = 
c 

controller  gain. 

11 


The  fractions  G1(s)/G2(s)  and  H^sJ/H^s)  are  equivalent  trans¬ 
fer  fractions  of  the  system  and  are  represented  by  ratios  of  polynomials 
which  upon  expansion  can  be  put  into  the  form 


n-k 


|i=Q 

m+n-i 


mfn-  i 
s 


> 


or  in  factored  form 


v  y 

sk  n  (s+z  )  n 
M>=1  ^q=v*l 


[«*(o  +>  )]  [st(o  )J 


s\!i(s+PiL!U  [o+(vj"t)]  [s+(°t->t)] 


where 


n-k-v  _  mfn-jfc  -g 
~~2~~  ’  r  =  - 2 - 


(1) 


(2) 


where  the  a's  and  b's  are  real,  the  poles  and  zeros  can  be  either  real 
or  complex  in  conjugate  pairs,  and  m  >  1. 

There  are  two  transfer  functions  of  fundamental  importance. 
OPEN  LOOP  TRANSFER  FUNCTION: 


CLOSED  LOOP  TRANSFER  FUNCTION: 


The  linear  system's  closed-loop  characteristic  equation, 


which  determines  stability  is 


KcKpG1(s)H1(s>  +  G2(s)H2(s)  =  0. 


(5) 


Dividing  both  sides  of  Equation  (5)  by  G„( s )H~( s)  gives 


Gj_(s)H1(s) 


KcKp  G*(s)H*(s)  +  1  *  0 


(6) 


If  we  let 


K  K  =  K* 
c  p 


Gi(s)Hi(s)  =  A(s),  and 


GgCsjH^s)  =  B(s),  then  Equation  (6)  becomes 


(7) 

(8) 

(9) 


(10) 


In  terms  of  the  system  open-loop  transfer  function, 


E 


+  1=0. 


(11) 


This  expression  relates  the  system's  closed-loop  characteristic  equa¬ 
tion  to  the  system's  open-loop  transfer  function. 

« 

The  steady  state  frequency  response  of  the  open- loop  transfer 
is  obtained  by  evaluating  Equation  (3)  along  the  imaginary  axis  for 
values  of  u)  between  zero  and  infinity  rad/s. 


The  result  is  a  vector  quantity, 


F 

E 


O  <  lb  <  oo 


K»  A(s) 

B(s7 


O  < 


S  <  oo 


(12) 


13 


where 


=  scalar  magnitude  of  open-loop  transfer  function  between 
0  <  it  <  ®  ,  and 

<h 

oL  =  direction  of  magnitude  of  open-loop  transfer  function 
between  0  <  iw  <  ®  . 

Equation  (10)  can  be  rewritten  as 


M  .  _  +  U0°  =  0, 

oL  ^  oL  * 

(13) 

or 

M  T  Z-0  _  =  1  ZLl80°. 
oL  oL 

(14) 

Equation  14  is  the  form  required  to  provide  quantitative  stability 
information  for  a  closed-loop  system. 

There  are  two  figures  of  merit  which  provide  a  measure  of 
system  stability.  The  only  restriction  is  that  all  roots  of  the  system's 
open- loop  transfer  function  denominator,  B(s),  be  <  0. 


GAIN  MARGIN  is  the  amount  that  M  T  must  be  increased  or  de- 
■  OL 

creased  to  make  it  equal  to  1,  only  when  the  orientation  of  magnitude 
of  the  open- loop  transfer  function,  is  l80°. 


PHASE  MARGIN  is  the  amount  that  0  T  must  be  increased  or  de- 

-  O  oL 


creased  to  make  it  equal  to  180  ,  only  when  the  magnitude  of  the  open- 


loop  transfer  function, 


These  concepts  can  be  related  directly  to  Equation  14. 
Incipient  instability  is  identified  with  either  a  zero  gain  or  zero 
phase  margin. 


The  scalai  magnitude  of  the  open- loop  transfer  function  is 
generally  expressed  in  decibels.  The  decibel  equivalent,  (db),  of 

a  number,  Mq^,  is  (db)  =  20  lo810  MoL* 
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The  gain  margin  and  phase  margin  of  a  linear  closed-loop 
system  as  shown  in  Figure  1  and  Equation  (4)  can  be  readily  determined 
by  analyzing  the  open- loop  transfer  functions  of  Equation  (3)  in  accord¬ 
ance  with  Equation  (l4). 

2.2  Nonlinear  System. 

V.M.  Popov's  plots  of  system  transfer  functions  give  some 
insight  into  stability  of  nonlinear  systems. Figure  2  shows 
the  system  configurations  considered  in  this  analysis.  It  contains  a 
linear  part  and  one  nonlinear  element  and  is  subject  to  the  following 
restrictions. 


pigure  2.  Nonlinear  Closed  Loop  System 


1.  The  roots  of  G^Cs)  have  negative  real  parts  and  simple  or 
multiple  roots  on  the  imaginary  axis  (one  root  pennitted  at  origin). 

2.  (V) 

3.  The  linear  system  of  Figure  2,  obtained  by  substituting 
0=co,  is  stable. 
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If  there  exists  a  non-negative  real  number  q  such  that 


[(1+jqt)  G(ju>)]  +  -|- 


>  0 


(16) 


for  all  u>>0,  where  G  is  defined  by 


Gi(s) 

S^sj’ 


the  nonlinear  system  is 


asympotically  stable  in  the  large. 

Since  G  has  been  defined  as  the  ratio  of  two  polynomials,  it 
may  be  expressed  as 


G(>)  = 


G1R(u>)  +  JG1I(») 

G2r(“)  +  ^ 


(17) 


where 


where 


G^R,  G^p  =  real  parts  of  G^,  and  G^,  respectively,  and 

^11*  ^21  =  Parts  of  G^  and  G^,  respectively. 

Redefining  G(juu),  we  have 
G(ju>)  =  Gr(u>)  +  jGjCuu), 


_  G1RG2R  +  G1I  G2I 

gR  =  (g2R)^  +  (g2I)4  >  and 


(18) 


becomes 


GI  = 


G1I  °2R  *  °1R  °2I  . 
(g2RJ2  +  { g2I;2 


As  a  result,  the  inequality, 


Re  [(1  +  jq»)  0  (Ji>)]  +  -i-  >  0  , 
Gr(»)  -  ^(o.)  +  -i-  >  0. 


(19) 


(20) 


1 6 


define 


To  aid  in  a  graphical  interpretation  of  these  results,  we 


U  =  and 

(21) 

W  =  itGj(uu). 

The  inequality  of  Equation  (20)  may  be  written  as 

U.qW  +  y>0, 

where  (22) 

q  >  0,  K  >  0  . 

If  we  replace  the  inequality  of  Equation  (22)  by  an  equal  sign, 

U-qW  +  y=  0.  (23) 


Then  it  may  be  seen  that  inequality  of  Equation  (22)  repre¬ 
sents  those  points  in  the  U-W  plane  of  the  modified  phase  amplitude 
characteristic  which  are  to  the  right  of  the  line  of  Equation  23,  as 
shown  in  Figure  3. 

MODIFIED  PHASE- 
AMPLITUDE 
CHARACTERISTIC 

- U 


Figure  3.  Stability  Boundary-Popov  Modified  Phase 
Amplitude  Characteristic. 

The  nonlinear  system  synthesis  problem  will  be  concerned  with 
determining  the  variation  of  K  for  a  given  linear  G(s).  Figure  4  shows 
the  boundary  constraints  of  the  nonlinearity  0  (o  )/a  of  the  system  shown 
in  Figure  2. 
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Figure  4.  System  Nonlinearity. 

To  maintain  symmetry  with  the  Bode  Plots,  the  modified  magni¬ 
tude  and  phase  characteristics  are  plotted  on  identical  axes.  Values 
of  the  modified  amplitude  associated  with  the  modified  phase  character¬ 
istic  for  90°  <  0  <  270°  (-90°)  establish  the  contour  described  in 
Figure  3. 

Once  this  contour  is  established  in  the  U-W  plane,  the  maximum 
K  obtainable  is  determined  by  inspection  from  the  intersection  of  the 
U  axis  and  a  straight  line  tangent  to  the  modified  phase  amplitude 
characteristic . 

3.  GRAPHICAL  METHODS 

An  on-line  printout  of  the  Bode  diagram  for  an  open-loop 
transfer  function,  F/E,  or  for  a  closed-loop  transfer  function,  C/R, 
is  generated  for  the  linear  system.  Similarly,  printouts  of  the  Popov 
"modified  amplitude -phase  characteristics"  are  generated.  The  phase 
and  amplitude  axes  are  scaled  in  degrees  and  decibels  with  an  arbitrar¬ 
ily  selectable  scale  between  +  0  max  degrees  +  M  max  db,  respectively. 
The  frequency  axes  are  logarithmically  scaled  in  radians  per  second  over 
a  range  of  6  decades,  beginning  with  an  arbitrarily  selected  lowest 


18 


decade.  Values  of  amplitude  and  phase  are  calculated  for  55  uniformly 

distributed  frequencies  over  the  6-decade  range,  plus  17  additional 

A(  s  ) 

points  centered  about  each  root  of  '  B(  ^)  *  This  addition?!  dis¬ 

tribution  of  printouts  expressed  as  a  multiple  of  each  frequency  is: 
(0.01,  0.05,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8, 

0.9,  1.0,  2.0,  4.0.  6.0,  8.0,  10). 

An  option  of  the  program  permits  the  density  of  the  Bode  plot 
to  be  either  sparse  or  dense.  The  sparse  plot  requires  K*  =  0  and  only 
one  set  of  frequencies  are  available,  from  which  the  amplitude  and 
phase  of  A(s)/B(s)  are  calculated.  The  dense  plot  permits  K*  to  have 
values  consistent  with  the  root  locus  subroutine. 

For  each  value  of  K*,  a  different  set  of  frequencies  will 
exist  for  K*  A(s)  +  B(s). 

A  larger  number  of  amplitude  and  phase  values  will  be  gener¬ 
ated.  Use  of  the  dense  plot  capability  may  be  necessary  to  obtain 
precise  definition  of  the  graphical  plot  in  some  cases. 

Any  number  of  simultaneous  expands  of  different  portions  of 
the  6-decade  frequency  range  may  be  plotted. 

4.  DATA  FORMATS 

4.1  General  Description. 

This  is  a  subroutine  designed  to  supplement  the  methods  pre¬ 
sented  in  AMSAA  TM  No.  21,  A  Linear  Closed  Loop  System  Analysis  Proce¬ 
dure  Using  Line  Printer  Plots  of  Characteristic  Equation  Root  Loci. 

This  subroutine  utilizes  the  polynomial  multiplication  and  root  locus 
methods  discussed  in  AMSAA  1M  No.  21,  and  extends  the  analysis  capa¬ 
bilities  to  include  the  following  available  options: 

1.  Root  Locus  Plots. 

2.  Frequency  Response  Plots. 

3.  Popov  Modified  Phase-Amplitude  Characteristics  Plots. 


19 


For  the  root  locus  option,  the  output  is  identical  with  the 
output  obtained  from  AMSAA  TM  No.  21. 

For  the  frequency- response  option,  the  output  is  a  log-log 
plot  cf  steady  state  amplitude  vs  frequency  and  a  semi- log  plot  of 
steady  state  phase  angle  vs  frequency.  Linear  expands  of  amplitude 
and  phase  for  any  frequency  range  may  be  plotted. 

For  the  Popov  modified  phase-amplitude  characteristic,  the 
output  is  a  log-log  plot  of  steady  state  amplitude  vs  frequency  and 
semi-log  plot  of  steady  state  phase  angle  vs  frequency.  Linear  ex¬ 
pands  of  amplitude  and  phase  for  any  frequency  range  may  be  plotted. 

4.2  Input  Description. 

a.  Program  option  control  card  that  precedes  all  data  (see 

below). 

b.  Program  option  control  card. 

To  run  the  program  in  its  entirety  with  a  dense  plot  in  the 
frequency  response,  i.e.:  Root  Locus  plus  Expands  and  Frequency 
Response  and  Popov  Plots  plus  expands  place  a  one  (l)  in  columns 
10,  20,  30,  and  40. 

To  run  only  the  root  locus  log  plot,  place  a  one  (l)  in 
column  10  and  zeros  (0)  in  columns  20,  30>  40. 

To  run  the  root  locus  expand  plot,  place  a  one  (l)  in  column 
20  and  zeros  (0)  in  columns  10,  30,  40. 

To  run  the  sparse  frequency  response  Popov  plot  only  place  a 
one  (l)  in  column  30  and  zeros  in  columns  10,  20,  and  40. 

To  run  the  dense  frequency  response  plot  only,  place  a  one  (l) 
in  columns  30  and  40  and  zeros  in  columns  10  and  20. 

To  void  all  plot  routines  place  a  zero  in  columns  10,  20,  30? 

and  40. 

c.  Comment  Card  that  follows  the  program  option  control  card 
for  identification  purposes. 
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d.  Program  Data  Cards. 

Same  format  as  described  in  Chapter  IV  (Data  Formats),  para¬ 
graph  A2  of  AMSAA  TM  No.  21.  Control  Card  used  to  call  All  Plot  Sub¬ 
routines  is  same  as  control  card  used  to  call  log  plot  subroutine  of 
AMSAA  TM  No.  21  (Section  IV,  paragraph  B2). 

Control  cards  used  to  call  root  locus  option  are  same  as  control 
cards  used  to  call  root  locus  expand  option,  paragraph  B2,  section  IV  of 
AMSAA  TM  No.  21.  Table  1  demonstrates  this  input  format  for  the  dense 
frequency  response  option.  Control  cards  used  to  call  frequency  response, 
Popov  plots,  and  expands  of  both  are  additional  features  of  this  report. 


DESCRIPTION 

COLUMNS 

DATA 

Data  appears  only  once 

a.  -  1-10 

Smallest  frequency, 
rad/s  (real) 

b.  -  11-20 

Max  phase  angle,  deg, 

(real) 

c.  -  21-30 

Max  amplitude,  DB  (real) 

Frequency  response,  Popov 

a.  -  1-10 

Zero, if  no  expands  desired 
(integer), 

One,  if  expands  desired 

Number  of  expands  desired 

a.  -  1-10 

No.  of  freq,  response 
phase  angle  expands  (inte¬ 
ger) 

b.  -  11-20 

No.  of  freq  response  ampli¬ 
tude  expands  (integer) 

c.  -  21-30 

No.  of  Popov  phase  angle 
expands  (integer) 

d.  -  31-40 

No.  of  Popov  amplitude 
expands  (integer) 

Region  of  expand,  one  card  for 

1-10 

Smallest  freq,  rad/s 

each  expand  in  consecutive 
order  with  a,b,c,d,  above. 

(real) 

11-20 

Largest  freq,  rad/s 
(real) 
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4.3.  Output  Description 


Depending  on  the  options  specified,  the  output  consists  of: 

Poot  Locus  Plot  Option 

a.  A  list  of  all  of  the  roots  plotted  on  the  complex 
frequency  plane  log  plot. 

b.  Two  log  plots  of  the  roots  calculated  by  the  polynomial 
multiplication  and  root  locus  program, 

c.  The  specified  number  of  linear  expand  plots  of  selectable 
regions  in  the  complex  frequency  plane. 

FREQUENCY  RESPONSE /POPOV  PLOT  OPTION 

a.  A  list  of  amplitude,  amplitude  in  decibels,  and  phase 
angle  in  degrees  for  various  frequencies  specified  in 
Section  3. 

b.  A  log- log  plot  of  amplitude  vs  frequency  and  a  semi-log 
plot  of  phase  angle  vs  frequency. 

c.  The  specified  number  of  linear  expands  of  the  frequency 
axis  for  the  amplitude  and  phase  angles  of  the  specified 
transfer  function. 


5.  EXAMPLES 

The  two  examples  discussed  in  AMSAA  TM  No.  21  will  be  analyzed, 
using  the  frequency  response/Popov  plot  option.  The  root  loci  of 
A(s)/B(s)  using  K*  as  a  parameter  have  previously  been  determined.  The 
characteristic  equation  of  the  system  shown  in  Figure  2  of  Reference  1 
is 

K*A( s )  +  B(s)  =  0. 
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The  open-loop  transfer  function  for  this  expression,  as  given 
in  Equation  (12)  is 


EXAMPLE  1. 
where 

K*  =  -0.005 

A^s)  =  2092 (S^+1.2S- 1672  )(S2+3.2S-*-25, 600 '+491S2(S2+6S+^  00) 

A2(s)  =  0.08l(464)(S+2.7)(S2+3.2S+25,600) 

A3(s)  =  0. 95(464)  (S+2. 7)  (S2+3.2S+25, 600) 

B(s)  =  S(S2+6S+400)(S2+3-2S+25,600) 

Table  1  shows  the  input  card  format  for  this  open- loop  trans¬ 
fer  function. 

Lines  2  through  24  apply  to  the  polynomial  multiplication  and 
root  locus  methods  and  lines  26  through  36  apply  to  the  frequency 
response  and  Popov  condition  plot  subroutines. 

The  basic  tabulated  numerical  output  is  shown  in  Tables  2  and 
3  and  is  similar  to  Tables  2  and  3  of  Reference  1.  A  tabulation  of  the 
open-loop  frequency  response  gain  and  phase  and  Popov  modified  phase 
amplitude  characteristics  are  shown  in  Table  4.  Figures  5  and  6  are 
the  Bode  plots  of  F/E  and  Figures  7  and  8  are  the  Popov  plots  of  F/E. 
Two  regions  are  expanded  to  demonstrate  capability  of  this  program. 
First  the  regions  of  frequencies  between  1.0  rad/s  and  300  rad/s 
are  shown  in  Figures  9  and  10  (frequency  response),  11  and  12  (Popov 
response)  and  secondly  the  region  of  frequencies  between  120  rad/s 
and  200  rad/s  shown  in  Figures  13  and  l4  (frequency  response),  15  and 
16  (Popov  response)  are  plotted.  Slope  information  consistent  with 
Bode  principles  may  be  applied  directly  to  these  plots. 
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TABLE  1.  INPUT  FOR  EXAMPLE  1 


COLUMNS 

1-10  , 

11-20 

21-30 

31-40  Remarks 

0 

0 

1 

1  Program  option 

control  card 

Frequency  Response 

Comment  l 

■  card 

0 

^  -9 

l.xlO  y 

1  -9 

l.xlO  y 

504  Problem 

data  cards 

0 

1  w  CO 

4 

3 

4 

H  *1  0 

<  e 

v#  H 

*5 

O  H 

rt> 

p  0 

►xj 

i 

.005 

P 

0 

3 

25600. 

3.2 

1  S3 

<g  3 

3 

-1672 . 

1.2 

1  P 

-L.  H 

^  p 
p  «+ 

1 

-2092 . 

s 

*3 

jr  p 

1 

.005 

p 

CO 

te> 

3 

0. 

0. 

-491.5  £ 

3 

• 

0 

0 

-3- 

6 • 

1 

-L. 

P 

1 

.081 

H* 

O 

£ 

2 

.  0. 

464. 

P 

c* 

2 

2.7 

1. 

H* 

O 

525 

3 

25600. 

3.2 

1.  3 

O 

• 

1 

.95 

§ 

ro 

1 

464. 

a 

H 

v« 

2 

2.7 

1. 

y 

0 

3 

25600. 

3.2 

1.  §• 

S' 

p 

2 

0. 

1. 

1  s 

-L  •  O 

c+ 

3 

• 

O 

O 

6. 

0) 

1  ^ 

3 

25600. 

3.2 

-L  •  W 

0 

0 

0 

10000  Same  as  paragraph  B2, 

Chapter  IV,  TM  No.  21, 
Control  card  used  to 
call  all  plot  sub- 

180. 

routines . 

.01 

100. 

CO  h(j 

1 

E5  2 

2 

2 

2 

2 

q  *0 

0  <*  c 

1. 

300. 

P  CO 

c+  P 

120. 

200. 

H*  P  O 

3  m 

1. 

300. 

P  *3 

•  O 

120. 

200. 

0  CD 

CO  CO 

1. 

300. 

n> 

0 

120. 

200. 

3 

H  M 

1. 

300. 

0  <0 
c+ 

120. 

200. 

1  P 

3 

r  Cb 
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EXAMPLE  2 


Table  6,  and  Figures  12  through  l4  of  Reference  1  are  the 
A(s)/B(s)  for  a  24th  order  linear  system.  Table  5  gives  the  F/E 
data  input  for  this  system,  and  Table  6  lists  open- loop  transfer  func¬ 
tion.  Figures  17  and  18  show  the  Bode  plot  and  Figures  19  and  20  show 
the  Popov  plot.  An  expand  plot  of  the  region  between  1  rad/s  <  U)  <  200 
rad/s  is  shown  in  Figures  21  and  22  for  the  frequency  response  and  in 
Figures  23  and  24  for  the  Popov  response. 

These  two  examples  demonstrate  the  significant  systems  analy¬ 
sis  capability  that  can  be  realized  by  using  the  tools  developed  in 
this  document  and  the  related  programs  contained  in  References  1  and  2. 
Detailed  conclusions  can  be  realized  when  high  order  feedback  control 
systems  are  encountered. 

6.  CONCLUSIONS 

The  method  described  in  this  memorandum  for  computing  linear 
system  root  loci,  steady  state  amplitude  and  phase,  and  the  nonlinear 
system  steady  state  modified  phase-amplitude  characteristic  has  the 
following  merits: 

1.  Accepts  input  in  the  form  of  sums  of  products  of  poly¬ 
nomials  • 

2.  Produces  a  log  plot  of  the  entire  complex  plane  in  two 
graphical  displays. 

3.  Produces  selectable  scaled  linear  plots  of  regions  in 
complex  frequency  plane. 

4.  Relates  the  closed-loop  gain  to  the  graphical  display. 

5.  Produces  a  magnitude  (db)  vs  log  frequency  steady  state 
response  plot. 

6.  Produces  a  phase  (deg)  vs  log  frequency  steady  state 
response  plot. 

7.  Produces  selectable  scaled  linear  plots  of  magnitude 
and  phase  of  the  steady  state  frequency  response. 
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TABLE  3.  EQUIVALENT  SYSTEM  OPEN  LOOP  POLYNOMIALS 

AND  HOOTS,  EXAMPLE  1. 
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TAbLE  L.  STAIXARD  rREQUE !JCY  REGPOIJLE  AKD  PCPOV  RE3P0NJE, 

EXAMPLE  1. 
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TABLE  Continued 
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FIGUhE  6.  EXAMPLE  1,  LINEAR  SYSTEM  BODE,  MAGNITUDE 
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JULY*  I9tB9  BRLESC  FORTRAN  2  ANO  A 
»SA046  J.ANDRESE  BLOG  390  X261I  • 


TT  A.. 51 
OOT  LOCUS 


15 

t  MAX? ( 10) MINS 

DO  10  J*l,100  1 

CALL  RTLOCS  2 

10  CONTINUE  3 

ENO  A 

5 

6 
7 
a 

SUBROUTINE  RTLOCS  9 

DIMENSION  SAVE II 100* 100) , SAVE 2 ( 1 00, ICO )  ,XA( 2C00) ,XB (2000 1  1G 

DIMENSION  X(100),IGA(  100),IGB(  100  )  ,  C(  100 )  ,D(  100)  ,  ANS(  IOC  )  •  11 

1SAVE(100),ERA( 100 ) ,  AS( 100  )  *  OS ( 100)  ,A( 100) ,B( 100 ) ,ROGTR f ICO ) •  12 

2R00T I  ( 100)  »  ATM  100 )«  CK(  100  )  ,CKS( IOC)  13 

INTEGER  COMEN T( 20 )« DENSE, ROTLCS,  f XPND, FREQSR, STAR, STAR 1,STAR2*  l A 

•  STAR!, ONE, ZERO,  DOT, VEE * AYE, CASH, BLANK, 0,  15 

♦  P0INT(  130, 100),  SP0T(130,50)  16 

DATA  ONE/ 1 HI/ , ZERO/ 1 HO/ »D0T/1H*/»AY£/ 1HA/, VEE/ 1HV/ • BLANK /1H  /,  17 

♦  DASH/1H-7 , Q/ 1HI/*STAR1/1H0/*STAR 2/1 H*/, STAR 3/1 H*/  la 

LOGICAL  PRNT  19 

EQUIVALENCE  ( SPOT  1 1, l ) ,P0 INT ( l , 1 )),( STAR  1, ZERO ),< STAR3,D0T ) ,  20 

•  ( JACKIE, KJ,  II, ( JOANN, KK,J»  21 

COMMON/ FREER/ A, B  22 

COMMON  /INFOA/  0, BLANK, DASH  23 

COMMON  / INF07/  ICOUNT  ?A 

COMMON  SAVE1,SAVE2, POINT  25 

COMMON  / 1 NF08/ ROTLCS «  EXPNC, FREQSR  26 

COMMON  / INF09/  DENSE  27 

READ  (5,105)  ROTLCS,  EXPND, FREQSR,0ENSE  28 

READ!  *",111  COMENT  29 

WR!iE(6,13)  COMENT  30 

! r( ROTLCS* N£«w*ANO* EXPND *NE*G* ANU*FREQSR*NE*C )  GO  TO  37  31 

I F ( ROTLCS*  EC*  0* AND* EXPND*  EG *0 )  WRITc(6,3C)  32 

IF(FRCQSR*EQ*0*AND*FXPND*EQ*0)  WRITE(6,31)  33 

IF ( FREQSR* E 0*0* AN D*R0TLCS*EQ*0)  WRITE|6,32)  34 

IF (ROTLCS* NE*0*ANO*EXPND*NE*0)  WRITE(6,33)  35 

IF (FREQSR*NF*0*AND*EXPND*NE*0)  WRITE(6,34)  36 

1F(FREQSR*NE*0*ANO*ROTLCS*NE*0)  WRITt(6,35)  37 

38  D0277  IOU*l,20CO  38 

XA(I0U)«0.0  39 

277  XB(nU)*0.0  AO 

DO  A01  KK* 1 , ICO  A 1 

DO  AOl  K J* 1 , l 30  A2 

A01  POINT (KJ,KK)* BLANK  A3 

00  A02  KK*  1 ,130  AA 

PO INT (KK ,5 l )*OASH  A5 

402  POI NT (KK,50 ) *DASH  A6 

DO  A03  KK* 1,100  A 7 

00  A03  KJ*3A, 125, 13  A8 

403  POINT(KJ,KK)*Q  49 

00  404  K  J* 34, 124  50 

00  404  KK*  7 ,  93,7  51 

IF(KK*EQ*49)  KK*5A  52 

404  POINT (KJ,KK)*OASH  53 

00  405  KK*7,  93,7  54 

IF ( KK* EQ*49 )  KK*58  55 

405  P0INT(25,KK)*0NE  56 

00  406  KK*26,28  57 

406  P0INT(KK,7)*ZER0  59  58 


DO  407  KK«26,27  59 

407  P0INT(KK,14 ) *ZER0  60 

POINT (26,21 )*Z6R0  61 

00  408  KK*36,65,7  62 

IF (KK.EQ.49)  KK*5H  63 

408  P0INT(24,KK)*U0T  64 

P0INT(25,42)*ZER0  65 

POINT(25,58)*ZERO  66 

P0INT(26,42)*0NE  67 

P01NT(26,58)*0NF  68 

PO I  NT  1 26»  79) *Z  ERO  69 

PO I  NT ( 26»  86 ) *Z  ERG  70 

P0INT(27,ft6)*ZER0  71 

00  409  KK*26,28  72 

409  P0INT(KK,93)*ZER0  73 

00  420  KK* 1(100  74 

IF ( KK.EQ.4 )  KK*9R  75 

420  POINT (26*KK)*C  76 

P'JINT(26,4)=VCE  77 

P0INT(26»97)*AYF.  78 

00  7661  JACK  I  E*l  * 100  79 

00  7681  JOANN* It  100  -  80 

SAVE1 ( JACK  I Et  JOANN ) *0*0  81 

7681  SAVE 2 ( JACK IE, JOANN)* 0.0  82 

NO* 1  83 

ICOUNT-O  8A 

MI *0  85 

JZO* 1  «6 

Kl*l  87 

111*1  88 

140  FORMAT  I 10X • 7HCELTAK*  , 1 10, 5X , 19HP0LY. AODEO  IN  A  I S 1* •  1 10 » 5X, 19HP0LV.  89 

•  ADDED  IN  R ( S ) * • 1 10»5X*  9HPR0B.NC. *» 1 10  I  90 

12  READ  (5,101)  N,IA, IB, IPROB  91 

IF( (P.NE.OI  GO  TO  10  92 

IF (  I A,NE*0 )  GO  TO  10  93 

I F (  N.NF.O)  GO  TO  10  94 

IF (  IPROH.EQ.IGOOO)  GO  TO  1016  95 

101  FORMAT  ( 7 1 lu )  96 

10  WRI TE(6, 1 )  97 

00  5  I  *  l, 10t  98 

A ( I  I  •  0.0  99 

5  8(1)  *  0.0  100 

F0RKI*0.0  IC1 

F0RK2  *  0.0  102 

F0RK3  *  0.0  103 

WRITE  (6  , 102 ) IPROB  104 

102  FORMAT  (1H1,  9X,40HP0LVN0MIAL  MULTIPLICATION  ANO  ROOT  LCC(JS,44X,U  105 
1HPR0DLEM  NO., 15//)  106 

IF  (N)  21,  20,  21  107 

20  RE A0(5, 103)V, UV,VT  108 

WRITE(6,140)N, 1A, IB, IPROB  109 

WRITE(6,141)V,0V,VT  110 

141  FORMAT ( 1H0, 5X , 10HK-IN! T IAL*F20*10,5X, 12H INCREMENT  K«,F20. 10, 5X, 12H  111 

•K-TERMINATE*, F20. 10)  112 

103  FORMAT  (4F1G.0)  113 

GO  TO  25  116 

21  REA0(5,IG4)(X( I), 1*1, N)  115 

WRI TE (6, 142 )  (I,X(I),I-1,N)  116 

104  FORMAT ( 7E10.0 )  117 

25  READ  (5  ,105)(IGA  (I),  I  •  1,  IA)  60  118 
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WRITE  (6,143)  (  I,  IGA(  I),  1*1,  U)  U? 

143  FORMATI 10X.24HNUMBER  OF  POLY.  IN  GROUP.  I  r  NUMERATOR* *110 1  120 

142  FORMAT ( 10X.2HX ( , I3,?H)*,F20«10 )  121 

READ  15  *  105) ( IGBI J) «  J*  1.  IB)  122 

WRI TE (6 . 144 )  (J,IGB< J),J*l,t6>  123 

144  FORMAT  I iOX . 24HNUMBER  CF  POLY.  IN  GROUP.  I3.16H  OF  DENOMINATOR* *  1 10  124 

*)  125 

105  FORMAT (7110)  126 

300  00  15  I  *  1.100  127 

15  SAVF ( I )  *  0.0  128 

ISOEG  *  0  129 

JJ  *  1  130 

200  READ  (5  ,106)MDEG,(  C(I).  1*1.  6)  131 

IFIMDEG.GT.6)  GO  TO  48  132 

WRI  TE  (6.  145)  (I. C(I)»  IM.MOEG)  133 

IF  (MDEG  -  6)  49,  49,  48  134 

48  READ  (5  ,104)(C(I),  I  *  7,  MDEG)  135 

WRI TE (6, 145 )  ( I . C ( I), 1*1, ROEG)  136 

106  FORMAT ( I 10.6E10.0 )  137 

49  IF  (IGA(JJ)-l)  50,  51,  5C  138 

51  DO  5fc  I*  l,  MDEG  139 

56  ANSI  I )  *  C(l)  140 

IADEG  *  MDEG  141 

IF  (ISDEG  -  MuEG)  52,  52,  53  142 

53  INDEG  *  ISDEG  143 

GO  TO  68  144 

52  INDEG  *  MDEG  145 

GO  TO  68  146 

50  READ  (5  ,  lOb)NOEG,  (0(1),  I  *1,  6)  147 

IFINDEG.GT.6)  GO  TO  54  14d 

WRITE (6, 145)  (1,0(1), I*l,NDEG)  149 

IF  INDEG  -  6)  55,  55,  54  150 

54  READ  (5  ,104)1011),  1*7,  NCEG )  151 

WRI TE (6, 145)  ( 1,0(1), I*1,NCEG)  152 

55  IADEG  *  NDEG  ♦  MDEG  -l  153 

CALL  POLMPY  I C.MOEG, D, NDEG, ANS  )  154 

IGA  (JJ)  *  IGA(JJ)  *1  155 

IF  (IGA(JJ)  -1)  65,  65,  64  156 

64  DO  60  I  *  1,  IADEG  157 

60  C(I)  *  ANSI  I )  158 

MOEG  * IADEG  159 

GO  TO  50  160 

65  IF  (ISDEG  -  IADEG)  66,66,67  161 

66  INDEG  *  IADEG  162 

GO  TC  68  163 

67  INDEG  «  ISOEG  164 

68  CALL  POLAOD  I  SAVE , ISOEG. ANS , I AOEG, ERA )  165 

145  FORMAT ( IOX, 2HC( , I3»2H)*»F?0*10)  166 

WRITE(6,1)  167 

1  FORMAT! IH1//I  168 

6800  IF  IFRAI INDEG) )6B03,  6802,  6803  169 

6802  INDEG  ■  INDEG  -  l  170 

IF  (INDEG)  6B01,  6R01,  6000  171 

6801  INDEG  «  l  172 

6803  JJ  «  JJ  ♦  l  173 

DO  70  I  *  l,  INDEG  174 

70  SAVE ( I )  *  ERA(I)  175 

ISDEG  *  INDEG  176 

IA  *  IA  -l  6l  177 

IF  (IA)  201,  201,  200  178 


201  IF  (FORKl)  202,  202,  203  179 

C  SAVE  NUMERATOR.  160 

202  00  220  I  *  1,  ISDEG  161 

220  AU)  *  SAVE  1 1 )  182 

IDA  *  ISOEG  163 

IA  «  IB  184 

FORKl  >  1,0  183 

00  230  1*1,  IA  186 

230  IGA ( I )  *  IGB(l)  187 

C  START  DENOMINATOR  186 

GO  TO  300  189 

C  SAVE  DENOM.  190 

203  00  240  I  *  1,  ISOEG  191 

240  B(I)  *  SAVE(I)  192 

I OB  *  ISOEG  193 

WRITE  (6  , 1C9 )  194 

109  FORMAT  ( 10X.41HC0EFF IC IENTS  ARE  GIVEN  IN  ASCENDING  ORDER/////)  193 

339  IF  ( A ( IDA) )  340,  341,  340  196 

341  IDA  =  IDA  -  1  197 

IF  1  IDA )  345,  345,  339  198 

345  WRITC  (6  ,120  199 

120  FORMAT  ( 1H0, 10X, 20HP0LVN0MI AL  A  IS  ZERO//)  200 

F0RK2  *  1.0  201 

GO  TO  410  202 

340  IF  (IOA  -  2)  346,  347,  335  203 

346  WRITE(6,1?1)A( 1)  204 

STAR=STAR1  205 

PRNT*. TRUE.  206 

121  FORMAT  ( 1H0,  1GX, 28HP0LYN0MI AL  A  IS  A  CONSTANT  *,  1PU16. 7// )  207 

GO  TO  410  208 

347  ROOT  *  -  All)  /  AI2)  209 

WRITE  (6  , 133 ) A( 1),  AI2)  210 

133  FORMAT  (10X,21HTHE  COEFFICIENTS  OF  A/IP2E20.7)  211 

WRITF  16  , 12? ) ROOT  212 

STAR*STAR1  213 

PRNT*. TRUE.  214 

ANUMP1*R00T  215 

ANUMH?*0.0  216 

CALL  PLCTF  R( SAVE  1 , S  AVE2, ANUMB1 ,ANUMB2,P0 1NT, XA, XB,MI *PRNT, STAR ,  I  1 1  217 

*,Nto  218 

122  FORMAT  I 1HC, 10X,23HR0CT  OF  POLYNOMIAL  A  I S, 1PIE16. 7//)  219 

GO  TO  410  22 0 

C  WRITE  POLYS  221 

335  ID1A  *  IDA  -1  222 

WRITE  (6,107) I01A, (All), 1*1, IDA)  223 

K  *  IDA  22 4 

DO  800  I  *  1, IDA  225 

ASH)  *  AIK)  226 

800  K  *  K-l  227 

IDP2AM0A  *2  228 

IU2A*  2  *ID1A  229 

CALL  KULLFR  I  AS, I D1A,ROOTR,ROOTI )  230 

DO  805  I  *  1, 101 A  231 

SAM*  100.  *  AMAX1I ABSIROOTKI I ) ) , ABSf ROOT I f III)  232 

IF  (SAM  ♦  ABSI ROOTRI I ) ).EQ.  SAM)  RCCTRI 1 )•  0.0  233 

IF  (SAM  ♦  ABSIROOTK  If  ).EQ.  SAM)  RCOTI  ( I )  *  0.0  234 

805  CONTINUE  235 

400  WRITF  (6,111)  (ROOTRI I ) , ROOT  I  (  I), I«l, 1D1A)  236 

CALL  ERCHEKI ROOT  1 , 1 0 1 A )  237 

PRNT*. TRUE.  52  238 


STAR* STAR l 
00  2  III-I.ID1A 
ANUMB1*R00TR( III  ) 

ANUMB2-ABS(R00TI( III)) 

CALL  PLOTERI SAVE 1,SAVE2, ANUM61«ANUM62, POINT, XA,X6«M| ,PRNT«STAR,I I I 

•  ,N0) 

.  2  CONTINUE 

410  IF(B('LD8)  I  4 1 1  r  412,  411 
412  IDB  «  106  -  1 

IF  C I OB  I  445,  445,  410 
445  WRITE  (6  ,123) 

123  FORMAT  ( 1H0, 10X, 2OHP0LYN0MIAL  6  IS  2ER0//I 
IF  (F0RK2U2, 450,12 

450  F0RK3  «  1.0 
GO  TO  698 

411  IF  ( I OB  -  2)  451,  452,  459 

451  WRITE  (6  ,124)6(108) 

STAR«STAR2 

PRNT-. FALSE, 

STAR-START 

124  FORMAT  (1H0, 10X, 28HP0LYN0MI AL  B  IS  A  CONSTANT  ■, IP1E14.7//) 

GO  TO  698 

452  ROOT  «  -BID  /  B(?) 

WRITE  (6  ,134)6(1),  B(2) 

134  FORMAT  (10X,21HTHE  COEFFICIENTS  OF  B/1P2E20.7) 

WRITE  (6  , 125 )R00T 

STAR-STAR2 

PRNT-.TRUE. 

ANUMPl-ROOT 

ANUMB2-0.0 

CALL  PLOTERI SAVE  1, SAVE2,ANUM6 1 ,ANUMB2, POINT, XA,X6, MI ,PRNT,STAR, 1 1 1 

♦  , NO) 

PRNT-. FALSE. 

STAR-START 

125  FORMAT  ( 1H0, 10X, 23HR00T  OF  POLYNOMIAL  B  IS, IPIEI6.7// ) 

GO  TO  698 

107  FORMAT  (10X,42HTHE  COEFFICIENTS  OF  POLYNOMIAL  A  (ORDER  •  I3,1H)/  ( 
UP6E20.7)  I 
499  I 0 1 B  *  IOB  -1 

WRITE  (6,108) I01B,(B( I),  1-1,106) 

106  FORMAT  (////10X,42HTHE  COEFFICIENTS  CF  POLYNOMIAL  B  (OROER  •  13, IM 
1)/  (IP6E20.7)) 

K  *  IOB 

DO  801  I  >  1, ICR 
BS(I)  •  B(K) 

B01  K  ■  K-l 

IGP2R*  I  OR  •  2 
I02B  ■  2  *  1018 

CALL  MULLER  (  BS, !OIB,ROOTR,ROOTI ) 

00  806  I  »  1,1016 

SAM  ^  100,  ♦  AMAX1 ( ARS( ROOTR( I )),ABS(ROOTI( I ) ) ) 

IF  (SAM  ♦  ABS(ROOTR( I)).60.  SAM)  ROOTRI I )■  0.0 
IF  ISAM  ♦  ABS  ( ROOT  III))*  EQ*  SAM)  RCOTKI)*  0*0 
606  CONTINUE 

500  WRITE  (6,112) (ROOTR( I ) , ROOT II I ), I*  1,1016) 

CALL  ERCHEKI ROOT  1 , 1016) 

STAR-STAR2 

PRNT-.TRUE. 
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ANUMB2*ABS(R00Tt(  HID  299 

CALL  PL07ERISAVE 1  ,SAVE2, ANUM61 ,ANUMB2, POINT, XA,X8« MI »PRNT , STAR, 1 1 1  300 

•»N0)  301 

3  CONTINUE  302 

PRNT*. FALSE.  303 

STAR*ST ART  304 

111  FORMAT  (IHO.llX, 14HTHE  ROOTS  OF  A/  (IPIE20.7,6H  ♦!  • 1P1E14. 7, 1PI  303 

1E20. 7V6H  ♦ I  tlPlEl4.7tlPlE20.7t6H  ♦!  ,IPIE14.7))  306 

112  FORMAT  ( IHO, 1 lXt 14HTHE  ROOTS  OF  B/  (IPIE20.7,6H  M  , 1PIE 14. 7 , 1P1  307 

IE20.7V6H  ♦  I  ,  IP1E14.7, 1PIE20.7.6H  ♦!  ,IP1E14,7))  308 

698  IF  (FCRK2) 12,699* 12  309 

699  IF  ( FORK  3 ) 12, 6991 , 12  310 

6991  WRITE  16  , 102IIPR3B  ill 

MSHEET  *  5  312 

C  START  K  CALCULATIONS  313 

IF  INI  702,702,533  314 

533  00  550  I*  1,  N  315 

00  541  J*  1,  I  DA  316 

541  ATKIJ)  *  XII)  *  A(J)  317 

C  COMPUTE  ROOTS  OF  K  *  A  ♦  6  318 

IOC*  MAXOI  IDA,  I  OB  I  319 

CALL  POLADD  I ATK, IDA, B, IDB,CK  I  320 

IOS  *  IOC  321 

554  IF  ICKI IOS) 1555,  557,  555  322 

557  IOS  *  IDS  -  l  323 

IF  1 1 OS )  658,358,  554  324 

558  WRITE  16  ,129)X|  I  )  3?5 

129  FORMAT  I IHO , 10X, 36HP0LYNCMI AL  K*A  ♦  B  IS  ZERO  FOR  K  *, IP1616. 7//)  326 

CO  TO  550  327 

555  IF  (IOS  -  2)  559,  560,  561  328 

559  WRITE  (6  , 1 30 )CK I IOS I ,  XII)  329 

130  FORMAT  I IHO, 10X, 35HP0LYN0MI AL  K*A  ♦  B  IS  A  CONSTANT  *  ,1P1E15.7,  10  330 

1H  FOR  K  *  ,IP1E14.7//|  331 

GO  TO  560  332 

560  ROOT  *  -CK Cl)  /  CM2 )  333 

WRITE  (6  ,131) ROOT ,  XII)  334 

131  FORMAT  I  IHO,  10X, 18HR00T  OF  K*A  ♦  B  «  , IP1E15. 7, 10M  FOR  K  *  ,1P1E1  335 

14.7//)  336 

GO  TO  550  3?  7 

561  K  *  IOS  338 

DO  803  J  *  1, IOS  339 

CKSIJI  *  CKIK)  340 

803  K  *  K  -  1  341 

IDIC  *  IDS  -  1  342 

I0P2C  *  IDS  *  2  343 

I02C  ■  2  •  I01C  344 

CALL  MULLER  I  CKS, I OIC, ROOTR,ROOT I  I  346 

DO  807  J  *  l, IDIC  34b 

SAM  *  100.  •  AMAX  II ABSIROOTRI J ) ) , ABSIROOTI I J) ) )  347 

IF  (SAM  ♦  ABSIROOTRI J))*EQ,SAM)  RCOTRIJ)  *  0.0  348 

IF  (SAM  ♦  ABSIROOTI!  J)  I.EQ.SAM)  RCCTKJ)  •  O.C  349 

807  CONTINUE  350 

WRITE  (6,808) I01C,X( I  ),(CK( J), J-1.10S)  351 

BOS  FORMAT  I///1CX,4BHTHE  COEFFICIENTS  OF  POLYNOMIAL  K«A  ♦  B  (ORDER  *  362 

1 13,7H)  K  *  1PIE16.7/I 1P6E20.7  ) )  353 

545  WRITE  (6,115) I ROOTRI J ) ,ROOT 1 1  J ), J*l, ID  1C )  354 

CALL  ERCHEKIROOT I •  I01C )  355 

CALL  SAVER ( ROOTR, ROOT  I, ID1C,SAVEI, SAVE2,JZ0,K1 )  356 

115  FORMAT  I  1H0,9X, 16HR00TS  OF  K*A  ♦  B/( 1P1E20. 7,6H  ♦  I  ,IPIE14.7, IP1  357 

1E20. 7,6H  ♦  I  , IP1E14.7, IP1E20.7,6H  ♦  I  ,IPIE14.7))  358 
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5452  MSHEET  «  HSHEET  -  I  359 

IF  (MSHEET)  546,  546*  550  360 

546  WRITE  (6  »102)IPR0B  361 

MSHEfcT  *  5  362 

550  CONTINUE  363 

GO  TO  12  364 

702  00  705  J  *  1,  IOA  365 

705  ATR(J)  *  V  •  A(J)  366 

C  COMPUTE  ROOTS  OF  K  •  A  ♦  I  367 

IOC*  MAXOIIDA,  | OR )  368 

CALL  POLADD  (  ATK, I OA, R, 100, CK )  369 

IDS  *  IOC  370 

754  IF  (CM  IOS )  )  755*  757,  755  371 

757  IOS  *  IOS  -  l  ,  372 

IF  (IOS)  758,  758,  754  373 

758  WRITE  (6  ,129)Y  374 

GO  TO  711  375 

755  IF  (IDS  -  2)  759,  760,  761  376 

759  WRITE  (6  , 139ICKI IDS),  Y  377 

GO  TO  711  378 

760  ROOT  =  -CMl)  /  CM? )  379 

WRITE  (6  ,131) ROOT ,  V  380 

GO  TO  711  381 

761  K  *  IDS  382 

00  ft 04  I  *  l, IDS  383 

CKSII)  *  CM  K  )  384 

804  K  *  K  -1  385 

ID1C  *  IOS  -  1  386 

IDP2C  *  IOS  *  2  387 

I02C  *  2  *  I01C  388 

CALL  MULLER  I  CKS, ID1C,ROCTR*ROOT! )  389 

00  809  I  *  1,I01C  390 

SAM  *  100.  *  AMAXll ABSIROOTRI I )),ABS( ROOT t ( I ) )  )  391 

IF  (SAM  ♦  ABSIROOTRI I) ).EQ.  SAM)  ROCTR(I)*  0,0  392 

IF  (SAM  ♦  ABSIROOTK I) J.EQ.  SAM)  RCOTMI)*  0.0  393 

809  CONTINUE  394 

WRITE  (6,808)  I01C,Y,  (CM  I  ),  1*1, IOS)  395 

MRtTE  (6,115)  (ROOTRI J),R00TI( J),J*l, I01CI  396 

CALL  ERCHEKIROOTI, I01C)  397 

CALL  SAVER! KOOTR, ROOT  t, 1 01C,SAVE l, SAVE2, JZO,K 1 )  398 

711  Y  «  V  ♦  OY  399 

IF  (Y  -  YT)  712,712,12  400 

712  MSHEET  «  MSHFfcT  -  l  401 

IF  (MSHEET)  713,  713,  70?  402 

713  WRITE  (6  , 10? ) IPROB  403 

MSHEET  *  5  404 

GO  TO  70?  405 

37  WRI TE (6 , 36 )  406 

GO  TO  38  407 

1016  CALL  PLO TER ( SAVE  1 ,SAVE2, ANUMR1,ANUMB2, POINT, XA, XB,M| ,PRNT,STAR • 1 1 1  408 

*,N0 )  409 

1018  IF (EXPN0*EQ«0)  GO  TO  1017  410 

REAO( 5, 22 )L  411 

22  FORMAT! 110)  412 

IFIL.EQ.0t  GO  TO  1017  413 

CALL  EXPAND! L,XA,XR)  414 

1017  I F ( FKEQSK, EQ*0 )  RETURN  415 

CALL  FREQRS(XA,XB)  416 

11  FORMAT! 20A4 )  417 

13  F0RMAT(///////////,20X,20A4)  65  418 


30  F0RMAT|////2OX,44H0NLY  THE  FREQUENCY  RESPONSE  HAS  BEEN  PLOTTEO  )  419 

31  FORMAT! ////20X, 45H0NLY  THE  ROOT  LOCUS  LOG  PLOT  HAS  BEEN  PLOTTEO  I  420 

32  FORMAT! ////20X,39H0NLV  THE  LINEAR  EXPANO  HAS  BEEN  PLOTTEO  I  421 

33  FORMAT ( ////20X,59HTHE  ROOT  LOCUS  LOG  PLOT  ANO  LINEAR  EXPANO  HAVE  B  .422 

•EEN  PLOTTEO  I  423 

34  FORMAT I ////20X,62HTHE  LINEAR  EXPANO  ANO  THE  FREQUENCY  RESPONSE  HAV  424 

•E  BEEN  PLOTTED  I  425 

35  FORMAT  I ////20Xf 6BHTHE  ROOT  LOCUS  LOG  PLOT  ANO  THE  FREQUENCY  RESPON  426 

•SE  HAVE  BEEN  PLOTTED  )  427 

36  FORMAT I ////23X,28HEVERYTHING  HAS  BEEN  PLOTTEO  )  428 

RETURN  429 

END  430 

C  431 

C  432 

C  433 

C  434 

SUBROUTINE  POLMPV  !A,N,B,M,C)  415 

DIMENSION  AU)vRII)9C!I)  436 

K  *  M4N  437 

DO  5  I *  I »  K  438 

5  C 1 1  I  *  0.0  439 

DO  10  I  *  l  *N  440 

L  *  1-1  441 

DO  10  J*  l »  M  442 

L  *  LM  443 

10  CIU  *  C(L  >  ♦  A!  I)*B|J)  444 

RETURN  445 

END  446 

C  447 

C  448 

C  449 

C  450 

SUBROUTINE  PGLAOO  (A,N,R,M,C)  451 

DIMFNSICN  A! I ) ,B( 1  ),CI 1 1  452 

IF  IN-M)  1*1,2  453 

1  NK  =  N  454 

GO  TO  5  455 

2  NK  =  M  456 

5  DO  10  1*1, NK  457 

10  C 1 1  I  *  A  (!) I  1 1  458 

NK  *  NK» 1  459 

IF  (N-M)  11,25,15  460 

11  00  20  I*NK,M  461 

20  C(I)  *  BID  462 

25  RETURN  463 

15  DO  30  I*NK,N  464 

30  C(I)  *  All)  465 

RETURN  466 

END  467 

C  468 

C  469 

C  470 

C  471 

SUBROUTINE  MULLER(C0E,N1,R00TR,R00TI)  472 

DIMENSION  COE! 1 ) ,R00TR(  1),R00TI( 1)  473 

N2*NU1  474 

N4*0  475 

l«NUl  476 

19  IFICOEI I ) )9,7,9  477 

7  N4-N4+1  478 
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ROOTR(N4)*0. 

ROOT I CN4|*0. 

I»!-l 

IHN4-N1I19, 37*19 
9  CONTINUE 

10  AXR«0.8 
AXI«0. 

1*1 

Nl«l 

ALPIR*AXR 

aipu*axi 

H»1 

G0T099 

11  BETIR*TEMR 
BET  1 I*TEPI 
AXR*0.85 
ALP2R=AXR 
ALP2 1* AX  I 
M«2 

C0T099 

12  BET2R«TEPR 
BET2 I*TEP1 
AXR*C.9 
ALP3R=AXR 
ALP3I *AX  I 


M»3 

G0T099 

13  BET3R-TEPR 
BET3I *TEPI 

14  TE1=ALP1R-ALP3R 
TE2=ALPl 1-ALP31 
TE5*AIP3R-ALP2R 
TE6*ALP3 1 “AL P2  I 
TEP*TE6*TE5»TE6*TE6 
TE3MTE1*TE54TE29TE6I/TEM 
TE4* ( TE2*T E5- TE14TE6 l/TEM 
TE7-TE3M. 

TE9»TE3*TE3-TE44TE4 

TE10*2«4TE3*TE4 

UE15*TE7*BET3R-TE4*BET3I 

0E16»TE7*BET3I4TE4*BET3R 

TEll*TE3*BET2R-TE4*BET2l4BETlR-0E15 

TE12*TE3*BET2I^TE49BET2R^BETII-0E16 

TE7*TE9-1. 

TEl*TE9*RET2R-TE 10*BET2I 
TE2*TE94BET2 I 4TE 104BET2R 
TE13»TE1-BET1R-TE7*BET3R^TE 10*BET3I 
TE14«TE2-BET1I-TE7*BET3I-TE10*8ET3R 
TE 15*DE 15*TE3-DE 16*TE4 
TE16*DE154TE4*DE16*TE3 

TEl*TE 13*TE 13-TE 14*TE14-4#4( TE 11 *TE 15-TE 12*TE 161 
TE2s2.4TEl34T6l4-4**<TE124TEl54TE114TEl6l 
TEM  *  SORT  I TEl*TEl4TE2*TE2) 

IF(TFl) 113,113,112 
113  TE4  *  SORT  1.5  *  U  EP  -  TED) 

TE3*.5*TE2/TE4 
GO  TO  111 

112  TE3  *  SORT  1.5  *  CTEP  ♦  TEH) 

IF(TE2) 110,200,200 


110  TE3*-TE3 
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200 

TE4-.5*T62/TE3 

539 

ill 

TE7«TEI3»TE3 

540 

TES«TE!4*TE4 

541 

TE9-TE13-TE3 

542 

TE10*TE14-TE4 

543 

TE1*2.*TE15 

544 

TE2*2.*TE16 

545 

IF ITE74TE7»TE8*TE8-TE9*TE9-TEI0*TE 101204, 204, 205 

546 

204 

TE7-TE9 

54  7 

TES«TE10 

548 

205 

TEH*TE7*TE7»TE8*TE8 

549 

TE3«ITE1*TE74TE2*TE8|/TEM 

550 

TE4*|TE2*TE7-TE14T€8I/TEH 

551 

AXR*ALP3R*TE3*TE5-TE4»TE6 

552 

AXI«ALP3I*TE3*TE6*TE4*TE5 

553 

ALP4R* AXR 

554 

ALP4l«AX| 

555 

H«4 

556 

GO  TO  99 

557 

15 

N6*l 

558 

38 

IF  (ABS  (HELL)  ♦  ABS  (BELL)  -  UE-20)  18,18,16 

559 

16 

TE7  *  ABS  ( ALP3R  -  AXR )  ♦  ABS  ( ALP3I  -  AXI ) 

560 

IF  (TE7  /  (ABS  (AXR)  ♦  ABS  (AXI))  >  WE- 7) 18, 18, 1 7 

561 

17 

N3«N3M 

562 

ALPIP*ALP2R 

563 

ALP1 l*ALP2 1 

564 

ALP2R*ALP3R 

565 

ALP2WALP3I 

566 

ALP3R*ALP4R 

567 

ALP3 1  * ALP4 1 

568 

BET IR*BET2R 

569 

BET1  l*BET2I 

570 

BET2R«BET3R 

571 

9ET2I«BET3! 

572 

BET3R*TEPR 

573 

BET3I*TEHI 

574 

1  F  (  *43-100 )  14,  18,  18 

575 

18 

N4«N4* 1 

5  76 

ROOTRI N4 )*ALP4R 

577 

ROOTI (N4)*ALP4I 

578 

N3>0 

579 

41 

IF ( N4-N1 )30, 37,37 

580 

37 

RETURN 

581 

30 

IF  (ABS  (R00TMN4))  -  l.E-5 )  10,  10,  31 

582 

31 

GO  TO ( 32 , 10  )  , L 

583 

32 

AXR*ALPIR 

584 

AXI *-ALP 1 I 

585 

ALP1 l*-ALPl 1 

586 

M«5 

587 

GO  TO  99 

588 

33 

BETIR»TEMR 

589 

BETl l*TEHI 

590 

AXR*ALP2R 

591 

AX I *-ALP2 1 

592 

ALP2 1 *-ALP2 1 

593 

M*6 

594 

GO  TO  99 

595 

34 

BET2R*TEMR 

596 

BET2I*TEMI 

597 

AXR>ALP3R 

598 

AXI*-ALP3I 

ALP3I*-ALP3I 

L«2 

M*3 

99  TEMR*C0EII) 

rem  *0.0 

001001*1, Nl 
TEI«TEMR*AXR-TEMI*AX  I 
TCMI*TEMI*AXR«TEMR*AXI 

100  TEMR*  T F 1  ♦COE  I !♦ 1 1 

HELL*TEMR 
BELL=TEMI 

42  IF ( N4 )102, 103,102 

102  D01G 1 1  *  1 «N4 
TEM1»AXR-R00TK( I ) 

TEM2*AX l-ROOT III ) 
TE1*TEM1*T€M1*TEM?*TEM2 
TE2*ITEMR*TEM1»TEMI*TEM2I/TE1 
TEMI *  I TEMI •T6M1-TEMR4TEM2 )/TE  l 

101  TEMR*TE2 

103  GO  T0I11,12«13,15,33,34),M 

END 


SUBROUTINE  EXPLOT  1 XARRAV * IC IMEN, Y ARRAY, KO(MEN,BL |Mt T.LET ) 
INTEGER  PO I NT  1 130* 100 1 *GRAPHLI  120,36), STAR 
DIMENSION  XARRAV ( IOIMEN ) * VARRA VI KOIMEN ) 

DIMENSION  ENCRM1 ( 120 1* ENCRM2I 58) 

DIMENSION  SAVE  11  100* 1001  * SAVE21 100, 100) 

COMMON  SAVE1*SAVE2, POINT 
EQUIVALENCE  I  POINT  1 1,  1),  GRAPH.!  1,  l)  I 
DATA  STAR/IHa/ 

WRITE  16*1) 

READ!*, 10)01,02 
D3*RL IMIT 
U4*-PL I M I T 

IF (D2.GT ,Dl )  CALL  SWITCHI C1,D2 ) 

I F ( D4.GT • 03 )  CALL  SW  ITCHI C3,04 ) 

ICl*03 

JDIMEN-KDIMEN 

CALL  RLANKRIGRAPHL, 1 2C,  SB, 120 ) 

DO  22  1* 1 , JC IMEN 
X*X ARRAV ( I  ) 

V* VARRAY ( I  ) 

IFf X. EQ.O.C. ANO.Y.EQ.O.O)  GC  TO  22 
CALL  LINAR  | X, L IM  ITX, 120, C 1 ,D2 .ENCRM1 ) 

CALL  LINAR  I V,LIM|TY, 58,03,04, ENCRM2) 
1FILIMITX.EC.0.0R.LIMITY.EG.0)  GO  TO  22 
LIMITX*121-LIMITX 
GRAPHLI L IMI TX,LIM(TY )*ST AR 
22  CONTINUE 

WRITE  16, 30)  I  Cl, LET, (GRAPHLI I, 1), 1*1, 120) 

IC1 *04 

WRITE  1 6,  34  Ml  GRAPHLI  I,J),  1*1,  120),J«2,57) 

WRITE  16, 30)  IC1, LET, (GRAPHLI 1, 58 ), 1*1, 120) 

WRI TE (6, 35 )  02,01 


1  FORMAT ( 1 H 1 ) 


10  FORMAT (6F10,0,20X) 
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30  F0RMAT(1X,I4,A4,1X,120A1I  669 

34  FORMAT  I 10X, 120A1 )  660 

33  FORMAT ( 3X,  l PIE 10*  3, 7HR AD/ SEC* 9 IX, IP IE 10*  3,  7HRA0/SEC I  661 

61  RETURN  662 

END  663 

664 
663 
666 
66  7 

SUBROUTINE  LINAR  I V, L|M|  TV.  101  MEN, 0IMNS1,DIMNS2, ENCRMT  I  668 

DIMENSION  ENCRMTI I0IMEN)  669 

DIMENS«0IKNS 1-DIMNS2  670 

LIMITV«0  671 

J* 1 01 MEN-l  672 

A»J  673 

0ELTA«0IMENS/A  674 

6NCRMT ( 1)*0IMNS1  673 

ENCRMT ( ! 01  MEN I *01 HNS 2  676 

UO  10  1*2, J  677 

10  ENCRMT f I )» ENCRMT I l-l I -DELTA  678 

00  11  1*1, J  679 

11  IF  IV,  LE.  FNCRMT I  I ) • AND, V,GE, ENCRMT  I  1*11  ICC  TO  20  680 

GO  TO  40  681 

20  LIM|TY«|  682 

40  RETURN  683 

END  684 

683 

686 

SUBROUTINE  SM|TCH(X,YI  687 

TEST*Y  688 

V»X  689 

X-TEST  690 

RETURN  691 

END  69 2 

693 

694 

SUBROUTINE  ERCMEKlX.il  693 

DIMENSION  XI 1001  696 

DATA  ERRLIM/0.1E-8/  697 

00  10  J*1,I  698 

10  IFIABSUU)  I.LT.ERRLIM)  X|j|*0.0  699 

RETURN  700 

END  701 

702 

703 

704 
703 

SUBROUTINE  EXPAND! ILOVE, X,Y  I  706 

OIMENSION  SAVEll 100, 100 1 , SAVE21 100, ICO  I , X21 20001, V21 2000 1 , XI2C00 1 ,  707 

•  Yf  2000 1 ,Xl 1 2000),Y11 2000 1 ,Z f 1231 ,21 f 87)  708 

INTEGER  SPOT! 130,30) , DASH, C, BLANK, PC  INTI  130* IOC I  709 

OATA  C/1HI/,  CASH/1H-/, BLANK/IH  /  710 

DATA  APLUS/0,0/,BPLUS/0«0/,AM|NUS/0»0/»BMINUS/0*0/  711 

COMMON  / I NF04/  O, BLANK, DASH  712 

COMMON  /PAYNE/AM, AP,BM,BP  713 

COMMON  /PJ/X1, VI  714 

COMMON  SAVE1.SAVE2, POINT  713 

EQUIVALENCE! AP, APLUS 1,1 AM,AMINUS I , i BM,8MINUS) * IBP,BPLUS)  716 

EQUIVALENCES  11)  ,Zt  1)  »,IY1(  1),Z1(1I)  717 

EQUIVALENCE  I  SPOT ( l, 1 1 ,P0 INTI  1 , 1 1 1 , 1 1 1, 1 2 1  718 

70 


uouu 


00  40  J*l» (LOVE 

719 

00  50  M* 1 « 2000 

720 

XlCM)-0.0 

721 

50 

V|(M|*0«0 

722 

CALL  BLANKRI SPOT, 130,50, 1261 

723 

CO  54  M* 1,2000 

724 

X2(M|*0»0 

725 

54 

V2fM|*0.0 

726 

RE A0I5, 10  )  OMEGA, ENCRMT , S IGMA,  OEL TA 

727 

OMEGA- ADS< OMEGA) 

728 

10 

F0RMATI4F10.0) 

729 

A-ENCRMT 

730 

PE RCNT*A*0# 01 

731 

O-OMEGA 

732 

B«B*PERCNT 

733 

aplus*b*gmega 

734 

O-AMINUS 

737 

11-1 

738 

AMINU$*CMEGA>B 

735 

C-APLUS 

736 

00  41  L*l,2000 

739 

IF(V(L).LT,0)G0  TO  41 

740 

IFIVILUGT.CIGO  TO  41 

741 

23 

Xl(tl)*X(U 

742 

VII  ll)-Y(L) 

743 

Il-IUl 

744 

41 

CONTINUE 

745 

A-OELTA 

746 

901 

PtRCNT- A*0*0l 

747 

B- SI CM A 

748 

903 

R*B«PERCNT 

749 

BPLUS*SIGMA*B 

750 

BMINUS-SIGMA-B 

751 

C-BPLUS 

752 

O-BMINUS 

753 

12*1 

754 

910 

DO  42  L*l,2003 

755 

IFIC.GT.0,0)  GO  TO  60 

756 

IFIXKU.GT.DI  GO  TO  42 

757 

IFCXKD.LT.C)  GO  TO  42 

758 

GO  TO  25 

759 

60 

IFfXllU.GT.C)  GO  TO  42 

760 

IFIXICLI.LT.D)  GO  TO  42 

761 

25 

X2II2)*X1(L) 

762 

V2I  I2)*Y1IU 

763 

12*124 l 

764 

42 

CONTINUE 

765 

CALL  SPL 1 T I X2 , V2 • SPOT , APLUS • AM INUS9 BPLUS 9 6M INUS 1 

766 

40 

CONTINUE 

767 

RETURN 

768 

END 

769 

770 

771 

772 

773 

SUBROUTINE  SPLIT  1 X, V, SPOT, APLUS, AMINUS,BPLUS9BMINUSI 

774 

DIMENSION  X 1 ( 2000 ) 9 Y 1 1 2000 1 

775 

DIMENSION  X(2GOO),Y(2000)9Z< 1251, Zll 87) 

776 

INTEGER  SPOT  1 130,50) 

777 

COMMON  /PJ/X1,V1  71 

778 

noon 


EQUIVALENCE IXI '1),Z(1)),(Y1(1),Z1(1M 

A*APLUS 

B*AM|NUS 

C*A-B 

0*8 PLUS 

E*BMINUS 

G*A8S(0> 

HsABS(E) 

F*G-H 

0ELTA*C/ 124.0 
DIFF*F/50.0 
00  11  J*l,124 
Z(J)*A 

11  A*A-0ELT  A 
Z(125)*B 

IF(D.LT.0.0ID1FF*-0IFF 
00  12  J* 1  , 49 
Zll JI*D 

12  0*D-DIFF 
Z 1 1 50 1 *E 

CALL  BRAKUPIZ,X,V,SPGT,Z1) 

RETURN 

END 


SUBROUT  I NE  BR  AKUPI VV , X , V , SPCT ,  XX  I 
LOGICAL  SKIP* SKIP 1 

INTEGER  SPOT ( 1 30*501 •  XPT.VPT.STAR 

DATA  STAR/  lf*X/ 

DIMENSION  XI 2000 I ,  Y( 2000 I  ,YV| 125 1 , XX ( 871 
COMMON  /PAYNE/ AM,AP,BM,BP 

EQUIVALENCE  CAP* APLUS )• ( AM, AMINUS ) ,  (  6M,6MINUS I , (8P.BPLUS) 
SKIP*. FALSE. 

SKI  PI*. FALSE. 

00  4*  J*1 . 2G00 
00  43  1*1,125 
IFIYIJI.NE.O.OI  GO  TO  TO 
IF (XI J) .NE.0.0 )  GO  TO  30 
L*  J*  1 
M*J*6 

00  60  N*L»M 

IFIN.GT.2000)  GO  TO  60 
IF (YIN). NE. 0.01  GO  TO  30 
IF(XIN).NE.O.O)  GO  TO  30 
60  CONTINUE 
GO  TO  40 

30  IF(SKIP)  GO  TO  20 

IF(Y(J).LT.YY(  I)  )  GO  TO  20 
YPT*  I 

SKIP*. TRUE. 

20  IFISKIPDGO  TO  4fl 
IF! I.GT.50 )  GO  TO  43 
IFIBPLUS.GT.O.O)  GO  TO  10 
IF ( X( JI.GT.XXI I ) )  GO  TO  43 
GO  TO  11 

10  if « x ( ji.lt.xxi  in  no  to  43 

11  IF (BPLUS.LT .0.0)  XPT«I 
IF (BPLUS.GE.O. 1  XPT*5l-I 
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c 

c 

c 

c 


c 

c 


SKI  PI*. TRUE* 

I FI.NOT.SK IP ICO  TO  «3 
50  1*125 

SPOT! VPT.XPTI*STAR 
SKIP*. FALSE. 

SKIPI*. FALSE. 

40  IFf SKIPICD  TO  50 
43  CONTINUE 
40  CALL  RITEITfSPOTI 
RETURN 
END 


SUBROUTINE  RITEITISPOTI 
COMMON  / PAYNE/ AM. AP«  BP#  OP 

EQUIVALENCE!  AP. APLUS I  •  f  AM  v  AHINUS I «  f  BP.  BHINUS I  •  ISP.BPLUSI 

INTEGER  SPOT! 130.501 

IFf BP.LT.BMI  GO  TO  40 

A*RP 

B*BM 

MRlTEI4.il 

1  FORMAT  1 1HI.40X.43HC0MPLEX  FREQUENCY  PLANE  .RIGHT  HA  NO  OUAORANT  I 
GO  TO  30 

40  A*BM 
B*BP 

MRITEf 4.21 

2  FORMAT !1H|.40X.43HC0MPLEX  FREQUENCY  PLANE.  LEFT  HANO  OU AC RANT  I 
30  MRITEI6.i2IAP.AM 

12  FORMAT f  114X.|5H< - J-0PEGA./1X.F8.2.  1I2X.FB.2I 

MRITEI6.14IB 

14  FORMAT I 124X.5HSIGMA./126X. 1H1./126X. IHV./12IX.FB.2I 
MRITEI6. ill  SPOT 

II  FORMAT! IX. I30AII 
MRITEI6.15IA 

15  FORMAT! 60X.27HL1NEAR  EXPANC  PLOT! RAO/SEC  I  .34X.F0.2t 
RETURN 

END 


9 

10 

50 


40 

30 

20 


SUBROUTINE  SAVERfROOTR.ROOTI.10lC.SAVEl.SAVE2.J20.Kll 

01  MENS  I  ON  S  AVE 1  f  100.  100 1 .  SAVE  2 1 1 90.  ICO  I  .ROOTRC 1001  .ROOTI 1 100 1 

IF fKl 1 30.9.10 

Kl*l 


I01*I01C 

I01C*IDIC*!K1-1I 

IFf  ItllC.GE.  100IG0  TO  30 

DO  40  I2AP*K1. I01C 

I Z A* I ZAP- IK l- II 

SAVE1 f JZG. IZAPl*ROOTRf IZA I 

SAVE2! JZC.IZAP|*ROOTIf IZAI 

CONTINUE 

Kl*l01C*l 

GO  TO  20 

JZ0«JZ0*1 

I01C*I01 

GO  TO  9 

RETURN 

END 
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SUBROUTINE  PLOTERf SAVEl, SAVE2, ANUMB1, ANUMB2, POINT, XA, XB , MI .PRNT,  903 

*STAR,NBB,NO)  904 

DIMENSION  SAVE II  100*  100) t SAVE2I 100* 100) , XAI2000I ,XB(2000)  905 

INTEGER  POINT  1 130*  100) , DASH, BLANK, STAR, ROTLCS, EXPND* FREQ SR, OENSE  906 

COMMON  / INF07/  ICOUNT  907 

COMMON  / I NFOH/ ROTLCS ,EXPN0,  FREQSR  908 

COMMON  / INF09/  DENSE  909 

LOGICAL  PRNT  910 

IF ( PRNT )  GO  TO  43  911 

DO  40  NAP* 1 t 100  912 

DO  40  NBP*1, 100  913 

IFISAVEllNABvNBRI  141,42,41  914 

42  IFISAVE2(NAB,N6B )  141,40,41  915 

41  ANUMB1*SAVE1(NAB,NBR )  916 

ANUMR2*SAVE2INAB,NBB»  917 

ANUMB2* ABS ( ANUMB2  )  918 

43  A2Z* ANUM62  919 

A2Z«AZZ«100COU,0  920 

NZZ*AZZ  921 

IFINZZ.EC.O)  GO  TO  45  922 

50  NBB*NPB*1  923 

45  CALL  EXCUT E ( ANUMP1 , ANUMB2, POINT, XA,XB, Ml ,N0, STAR )  924 

IF  I DENSE,NE, 1 )  GO  TO  10  925 

ICOUNT* ICOUNT ♦ l  926 

GO  TO  11  927 

10  IF  C  PRNT )  ICOUNT* ICOUNT* 1  928 

11  IFIPRNT)  RETURN  929 

40  CONTINUE  930 

IFIRCTLCS.EQ.G)  RETURN  931 

70  CALL  NRITIT (XA,XB  )  932 

CALL  PREPARf POINT  )  933 

RETURN  934 

END  __  935 

936 

937 

938 

939 

SUBROUTINE  EXCUTF I ANUMB1, ANUMB2, POINT, XA, XB,MI ,N0, STAR  I  940 

LOGICAL  SKlPl »SK IP2  ,LESS  941 

INTEGER  ROTLCS, EX  PND, FREQSR  942 

DIMENSION  XA( 200C ) ,XB( 2000)  943 

INTEGER  POINT!  130,  100) , DASH, BLANK, STAR  944 

COMMON  / I NFG8/ ROTLCS, EXPNO, FREQSR  945 

DATA  K1/10/,K2/ 100/, K3/ lOOO/,K4/lGOOO/,NEGONE/-l/  946 

LESS**FALSE*  947 

1*0  948 

J*0  949 

XA(N0)*ANUMB1  950 

XB(N0)*ABS( ANUMB2  I  951 

NO-NOM  952 

IFIROTLCS.EQ.O)  RETURN  953 

IFCABSIANUM81 I ,GT. 10000,0)  GO  TO  50  954 

IF  I ABSI ANUMB2 ) ,GT« 10000,0 )  GO  TO  50  955 

IF  I ABS I ANUMB1 ) • EQ,0,0 )  GO  TO  500  956 

IF ( ABSI ANUMB2 ) •EQ,0,0)  GO  TO  500  957 

IFIABSI  ANUMB2WLT. 0.001)  GO  TO  50  958 

7h 


o  o  o  o 


IF ( ABS( ANUMBl )*LT* 0*001)  GO  TO  50  959 

500  CALL  SCALE! IKl«K2»K3»K4t NEGONE* I» ANUMB1»  ICONS* LESS* SKI FI )  960 

CALL  SCALE2(KltK2*K3tK4» J« JCONS* ANUM82»SKlP2)  961 

CALL  WPOINT ( J» JCONStL* ANUM62*SKIP2)  962 

CALL  SPOINT 1 I* ANUMBl* LESS*  ICONS*L*POINT»SKIPl* N0*XA#XB*STARI  963 

50  RETURN  964 

END  965 

966 

967 

968 

969 

SUBROUTINE  SCALE l C Kl ,K2* K3, K4, NEGONE, I . ANUMBl , ICONS .LESS , SKI P 1 I  970 

LOGICAL  LESS  971 

LOGICAL  SKIP1  972 

SKlPl**FALSE.  973 

ICONS*!  974 

1*1  975 

AKEEP* ANUMBl  976 

21  NUMBI*ANUMB1  977 

NUMB* I ABS ( NUMB 1 )  978 

IFf NUMB* E0*0) GOTO  12  979 

I F ( NUMB1 ) 31 .40.40  980 

40  I F ( NUMB1*GE* 10 ) GO TO  11  981 

GO  TO  50  982 

31  LESS**  TRUE*  983 

XXX*-NUMB1  984 

NUMBl*XXX  985 

GOTO  40  986 

12  GOTO! 1* 2*3*4) « I  987 

1  ANUMB1*AKEEP  988 

ICONS*Kl  989 

RK1*KI  990 

ANUMBl* ANUMBl *RK 1  991 

1*2  992 

GOTO  21  993 

2  ANUMB1*AKEEP  994 

I C0NS*K2  995 

RK2*K2  996 

ANUMBl* ANUMBl ARK 2  997 

1*3  998 

GOTO  21  999 

3  ANUMBl* AKEEP  1000 

IC0NS*K3  1001 

RK3*K3  1002 

ANUMBl* ANUMBl *RK 3  1003 

1*4  1004 

GOTO  21  1005 

4  ANUMB1*AKEEP  1006 

IC0NS*K4  1007 

RK4*K4  1008 

ANUMB1*ANUMB1ARK4  1009 

NUMB1*ANUMB1  1010 

I F ( NUMB1*  EQ*0 )  GO  TO  51  1011 

GO  TO  50  1012 

11  SKIP1**TRUE*  1013 

GO  TO  (6*7»8»9)»I  1014 

6  ANUMBl*AKEEP  1015 

IC0NS*K1  1016 

RK1*K 1  1017 

ANUMB1*ANUMB1/RK1  7<-  1018 


noon 


(•2  1019 

GOTO  21  1020 

7  ANUMBl*AKEEP  1  1021 

IC0NS*K2  1022 

RK2«K2  1023 

ANUMB1*ANUMB1/RK2  102* 

1*3  1025 

GOTO  21  1026 

8  ANUHB1* AKEEP  1027 

IC0NS*K3  1028 

RK3=K3  1029 

ANUMBI*ANUMB1/RK3  1030 

1*4  1031 

GOTO  21  1032 

9  ANUMB1*AKEEP  1033 

IC0NS*K4  1034 

RK4*K4  1035 

ANUMB1*ANUMB1/RK4  1036 

GO  TO  50  1037 

51  1*5  1030 

50  RETURN  1039 

ENO  1040 

1041 

1042 

1043 

1044 

SUBROUTINE  SCALE 2  (K1  »K2»K3«K4»  J»  JCONS»  ANllHB2» SK IP2 )  1045 

LOGICAL  SKIP2  1046 

SKIP2*. FALSE,  1047 

JC0NS*1  1048 

J*l  1049 

ANUMR2*ABS ( ANUMB2  )  1050 

BKEEP*ANUHB2  1051 

20  NUHB2*ANUMB2  1052 

IF ( NUMB2, EQ,Q ) GO  TO  10  1053 

I F ( NUKB2,GE, 10 )G0  TO  11  1054 

GO  TO  50  1055 

10  GOTO (It2f3*4) » J  1056 

1  ANUHB2*BKEEP  1057 

JC0NS*K1  1058 

RK1*K1  1059 

ANUMB2*ANUMB2*RK1  1060 

J *2  1061 

GO  TO  20  1062 

2  ANUMB2*BKEEP  1063 

JC0NS*K2  1064 

RK2*K2  1065 

ANUMB2* ANUMB24RK2  1066 

J*3  1067 

GO  TO  20  1068 

3  ANUMR2*BKEEP  1069 

JC0NS*K3  1070 

RK3*K3  1071 

ANUMB?*ANUMB2*RK3  1072 

J*4  1C73 

GO  TO  20  1074 

4  ANUMB2*BKEEP  1075 

JC0NS*K4  1076 

RK4*K4  1077 

ANUMB2* ANUHB2*RK4  76  1078 


I 


noon 


NUMB2*ANUMB2 
IFINUMB2.EQ.0)  GO  TO  51 
GO  TO  50 
11  SKIP2*.TRUE. 

GO  TO  (6f7f8f9),J 

6  JC0NS«K1 
ANUMB2»BKEEP 
RKUKl 

ANUMB2* ANUMB2/RK  1 
/  J*2 

GO  TO  20 

7  JC0NS=K2 
ANUHB2-BKFEP 
RK2*K2 

ANUMB7* ANUMB2/RK2 
J*3 

GO  TO  20 

8  JC0NS=K3 
ANUHB2*BKEEP 
RK3*K3 

ANUMB2* ANUMB2/RK  3 
J*4 

GO  TO  20 

9  JC0NS*K4 
ANUMB2-BKEEP 
RK4*K4 

ANUMB2* ANUMB2/RK4 
GO  TO  50 
51  J*5 
50  RETURN 
END 


SUBROUTINE  WPOINTt Jt JC0NS*LtANUMB2»SKIP2 > 
LOGICAL  SKI P2 


L*0 

IDELTA*0 

IF! ANUMB2.GE.9.2) 
I F ( ANUMB2.GE.8.0) 
IF  I ANUMB2«GE«6*9 ) 
IF ( ANUMB2*GE*  5*9 ) 
IF| ANUMB2.GE.5.0I 
IFIANUMB2*GE.4.2) 
I F  4 ANUMB2* GE« 3*5  ) 
I F ( ANUMB2«GE*  2*9 ) 
I F  4 ANUMB2«GE*2*4  ) 
IFUNUMB2.GE.1.9) 
I F  4 ANUMB2.GE. 1.5) 
IFIANUMB2.GE.U2) 
GO  TO  211 

200  IF  I ANUMB2*  LT« 9# 2  I 

201  IF4 ANUM82.LT.B.0) 

202  IFI ANUHB2*LT*6«9) 

203  I F ( ANUHB2*  LT«  5*9  ) 

204  IFI ANUMB2*  LT«5*0  ) 

205  IFI ANUMB2*LT*4*2 ) 

206  IFI ANUMB2.LT. 3*5  I 

207  I F  I  ANIIMB2*  LT«  2*9 ) 


IOELTA* 12 
GO  TO  200 
GO  TO  201 
GO  TO  202 
GO  TO  203 
GO  TO  204 
GO  TO  205 
GO  TO  206 
GO  TO  207 
GO  TO  208 
GO  TO  209 
GO  TO  210 

IDELTA* 11 

IDELTAMO 

IDELTA*9 

I0ELTA*8 

IOELTA-7 

I0ELTA«6 

I0ELTA«5 

IDELTA>4 


77 


n  o  o  o 


208 

IFI ANUMR2.LT.2.4)  I0ELTA*3 

1139 

209 

1FIANUMB2.LT.1.9)  IDELTA*2 

1140 

210 

IFI ANUMB2.LT.  1.5)  IOELTA* 1 

1141 

211 

IFIANUMB2.LT. 1*2) I DELTA*0 

1142 

I DELTA* IJ- I DELTA 

1143 

IFIJ.EQ.5)  GO  TO  50 

1144 

IF I SK IP? )  GO  TO  41 

1145 

IFI JCCNS.EQ. 1 ) GOTO  1 

1146 

IFI JCONS. EQ. 10 )GOTO  10 

1147 

IFI JCONS.EQ. 100) GOTO  100 

1148 

IF  I JCCNS.EQ. 1000) GO^O  1000 

1149 

IFI  JCONS.EQ.  lOOOOGOrO  100C0 

1150 

1 

L  =  73MDELTA 

1151 

GO  TO  40 

1152 

10 

L=86+ IDELTA 

1153 

GO  TO  40 

1154 

100 

L=99+ IDELTA 

1155 

GO  TO  40 

1156 

1000 

L* 1 12* I  DEL  T  A 

1157 

GO  TO  40 

1158 

10000 

L=126 

1159 

GO  TO  40 

1160 

41 

KCONS=JCONS*l 

1161 

ski P2*. false. 

1162 

IFIKC0NS.E0.1DG0T0  11 

1163 

IFIKCONS.Eti.lODGOTO  101 

1164 

IFIKC0NS.EQ.1G01 )GOTO  1001 

1165 

IFI  KCONS. EQ. 10001 ) GOTO  10C01 

1166 

11 

L*60*IDELTA 

1167 

GO  TO  40 

1168 

101 

L=47+IDELTA 

1169 

GO  TO  40 

1170 

1001 

L*34«IDELTA 

1171 

GO  TO  40 

1172 

50 

L*  125 

1173 

GO  TO  40 

1174 

1000 1 

L*34 

1175 

40 

RETURN 

1176 

END 

1177 

1178 

1179 

1180 
1181 


SUBROUT  I NE  SPGINTI I,  ANUMB  1,  LESS, ICONS, L.PCINT* SKIP1 , NO, XA,XB , STAR)  1182 
DIMENSION  X A ( 2C00 )  ,X6( 200C )  1183 
INTEGER  POINT! 1 30, 100) * DASH .BLANK, STAR ,0  1184 
INTEGER  STAR1*STAR2  1185 
DATA  STARI/1H0/.STAR2/1H*/  1186 
DATA  FLANK/1H  /*DASH/1H-/,  Q/lHl/  1187 
LOGICAL  SKIP1  .LESS  1188 
AKEEP*ANUMBi  1189 
ANUMB 1*ARS I AKEEP )  1190 
KOOL-O  1191 
IF I I.FQ.5)G0  TO  1000  1192 
IFI ANUMB l.GE. 7.1  )K00L*7  1193 
IF ( ANUMB l.GE. 5.0 )G0  TO  1011  1194 
IF  I ANUMB l.GE. 4.0 )G0  TO  1111  1195 
IFIANUMB1.GE.3.1)  GO  TO  1211  1196 
IFI ANUMBl. GE.2.3 )G0T013l 1  1197 
IF( ANUMB l.GE.  1*6 )  GO  TO  1411  7p  1198 


IF ( ANUMB1.GE. 1.0 ) GO  TO  1511 

1199 

GO  TO  4011 

1200 

1011 

IFIANUMBl.LT. 7.1 )K00L«6 

1201 

1111 

I F ( ANUMB l. L  T. 5.0 ) K00L*5 

1202 

1211 

IF (ANUMB1.LT. 4.0 )K00L»4 

1203 

1311 

IF  IANUMB1.LT.  3.  DKOOL-3 

1204 

1411 

IFIANUMBl.LT. 2.3 IKOOL-2 

1205 

1511 

IFIANUMBl.LT. 1.6 ) KOOL3 1 

1206 

4011 

NCRMNT«B-KOCL 

1207 

IF( LESS  I GO  TO  40 

1208 

IFfSKlPlIGO  TO  41 

1209 

IF! ICONS. EO. 1 ) GOTO  1 

1210 

IF ( I  CONS. EQ. 10 ) GOTO  10 

1211 

IF  1  ICONS. EQ.IOOIGOTO  100 

1212 

I F (  ICONS. EO.IGOOGOTO  1000 

1213 

1 

K00L  =  7UK00L 

1214 

GO  TC  50 

1215 

10 

K00L=64*K00L 

1216 

GO  TO  50 

1217 

100 

K00L=574K00L 

1216 

GO  TO  50 

1219 

1000 

K00L=50»K0OL 

1220 

GO  TO  50 

1221 

41 

LC0NS=IC0NS+1 

1222 

SKI  PI*. FALSE. 

1223 

IFILCCNS.EO.IDGOTO  11 

1224 

IFILCONS.EO.IODGOTO  101 

1225 

IF ( LCONS. EO. 1001 ) GOTO  1001 

1226 

11 

KOOL*784KOOL 

122  7 

GO  TO  50 

1228 

101 

KOOL-85*KOOL 

1229 

GO  TC  50 

1230 

1001 

K00L=92 ♦KOOL 

1231 

50 

LESS3. FALSE. 

1232 

IF ( PO INT 1 L* KOOL ) • EC. DASH )  GC  TO  51 

1233 

I F ( POI NT ( L« KOCL ) • EQ. G  )  GO  TO  51 

1234 

14 

I F ( PO 1 NT ( L  *  KOOL ) • EQ. ST  AR 1 )  GO  TO  15 

1235 

IF(POINT(L*KOOLI.6Q.STAR2)  GO  TO  15 

1236 

IF(STAR.EO.STARl)  GO  TO  51 

1237 

IFISTAK.EQ.STAR?)  GO  TO  51 

1236 

1F( POINT ILt KOOL ) «NE. BLANK )  GO  TO  70 

1239 

51 

POINT ( L.KOOL ) 3 STAR 

1240 

ANUMB1* AKEEP 

1241 

GO  TO  70 

1242 

15 

L3L*1 

1243 

IFIL.EO. 130)  GO  TO  51 

1244 

GO  TC  14 

1245 

40 

IF(SKIPI)  GO  TO  42 

1246 

IZ3 ICONS+2 

1247 

IF  I IZ.EQ. 3 ) GOTO  3 

1246 

IFIIZ.EQ.12IG0T0  12 

1249 

IF( IZ.EQ. 102) GOTO  102 

1250 

IF( IZ.EQ. 1002 (GOTO  1002 

1251 

3 

K00L=21+NCRMNT 

1252 

GO  TO  50 

1253 

12 

K00L328*NCRMNT 

1254 

GO  TO  50 

1255 

102 

K00L=35*NCRMNT 

1256 

GO  TO  50 

1257 

1002 

K00L=42*NCRMNT  79 

1258 

GO  TO  50  1259 

42  I Y*  ICONS+3  1260 

SK I Pl=*  FALSE*  1261 

IF1 IY*EQ*13)GOTO  13  1262 

IFI IY.E0*103)G0T0  103  1263 

IFC IY.EQ.1003IG0T0  1003  1264 

IF( I Y*EQ* 10003 1 GO TO  10003  1265 

13  K00L=14+NCRMNT  1266 

GO  TO  50  1267 

103  K00L-7+NCRMNT  1268 

GO  TO  50  1269 

1003  KOOL-NCRMNT  1270 

GO  TO  50  1271 

10003  PO I  NT  I L» 1 ) =ST  AR  1272 

GO  TC  70  1273 

70  RETURN  1274 

END  1275 

C  1276 

C  1277 

C  1278 

C  1279 

SUBROUTINE  PREPAR C PC INT I  1280 

INTEGER  POINT! 130,100), PRT  1281 

DATA  PRT  / 6/  1282 

WR I TE (6f 1 )  1283 

1  FORMAT! IHl)  1284 

WR I TE 1 6, 1 1 )  1285 

11  FORMAT!///)  1286 

WRITE! PRT • 14 )  1287 

14  FORMAT!  75X*8H|_0G  PL0T,/58X,42HC0MPLEX  FREQUENCY  PLANE, LEFT  HAND  CU  1288 

•  AO RANT  ,/ 75X.9H ( KAD/S  EC )  )  1289 

WR I TE (6« 12 )  1290 

12  FORMAT! 32X, 5H10000 ,9X,4H 1CC0, 9X, 3H1C0, 10X, 2H10, 12X, IHl , l IX ,2H. 1 . 10  1291 

*X,l4H.0l< - J-OMEGA* /  18X,11HMINUS  SIGMA)  1292 

13  FORMAT !32X,5H10000,9X,4H1000.9X,3H1CO,10X,2HIO,12X,IHI,11X,2H#1* 10  1293 

*X»  14H«01< - J-OMEGA*/  18X*UH  PLUS  SIGMA)  1294 

DO  50  1=1,130  1295 

50  POINT! I,50)=P0INT! 1.51)  1296 

WR ITE! PR T.10) POINT  1297 

10  FORMAT! IX, 130A1)  1298 

WR I TE !6, 1 3 )  1299 

WR ITE! PRT , 15 )  1300 

15  FORMAT! 75X,8HL0G  PL0T,/58X,43HC0MPLEX  FREQUENCY  PLANE,RIGHT  HAND  Q  1301 

•U ADR ANT ,/ 75X , 9H! RAD/ SEC )  )  1302 

RETURN  1303 

END  1304 

C  1305 

C  1306 

C  1306 

C  130  7 

C  1308 

SUBROUTINE  WRI TI T ( XA, XB )  1309 

DIMENSION  XA! 2000 )«X8(2000)  1310 

WRI TC !6, 1 )  1311 

1  FORMAT! IHl)  1312 

WRITE(6,U)  1313 

11  FORMAT !5X,48HTHE  FOLLOWING  ROOTS  ARE  PLOTTEO  ON  THE  LOG  PLOT,/  1314 

1 5X , 99HR00TS  AT  THE  ORIGIN  ARE  NOT  PRINTED  OR  PLOTTEO,  ROCTS  ON  THE  1315 

*  J-OMEGA  AXIS  ARE  NOT  PLOTTED*  ,  1316 

2»//16X,5HS!GMA,25X, 7HJ-0MEGA, // )  fto  1317 


00  77  120=1,2000  1310 

XC=XB(IZ0)  1319 

IF(XC.GT.O.O)  XC—XC  1320 

IF(XC.NE.O.G)  GO  TO  22  1321 

20  IF1XA! IZO) )22«40,22  1322 

40  K1*IZ0  1323 

K2*IZ0+12  1324 

DO  50  K*Kl ,K2  1325 

IFIXAIKI.NE.O.O)  GO  TO  77  1326 

50  IF!XB(K).NE.C.O)  GO  TO  77  132/ 

GO  TO  60  1328 

22  WRI TE !6, 10 ) XA( IZO ) «XC  1329 

10  F0RMAT(5X,F20.9, 10X,F20.9)  1330 

77  CONTINUE  1331 

60  RETURN  1332 

FNO  1333 

C  1334 

C  1335 

SUBROUTINE  FREQRS ( SIGH AX, CMEGAY)  1336 

DIMENSION  SAVE  1 ( 100» 100 1  *  SAVE 2 1 100,100), ACOEFF (lOO),BCOEFF(lOO),  133  7 

*  AAAAAA(99),BBB8BB!99),SIGMAX!2C00) » OMEGA Y I 2000 ) »OEC (9) *  1338 

*  OMEGAS (55), OMEGA! 17) ,W LOG ( 3000 ), PHIL  IN! 3000) »MCBLINC30CO  1339 

*  ) »PHIPLN(3000) »  MPCBLNI 3GCC ) , AL IM( 7 ) ,B IGSEM! 300C ) •  1340 

*  EMSUBS ( 3000 )  1341 

REAL  MDBLIN.MPCRLN  1342 

INTEGER  FRQEXP, I EXP! 4 1  1343 

INTEGER  POINT! 130* 10G)  1344 

COMMON/ FRFEK/ACOEF Ft  BCOEFF  1345 

COMMON/ INFG2/ALIM  1346 

COMMON/ INF03/0MEGA  134/ 

COMMON/ 1NF06/ EMSUBS  1348 

COMMON  / INF07/  ICOUNT  1349 

COMMON  SAVE it  SAVE2t POINT  1350 

EQUIVALENCE  ! WLOGil) .SAVE1!  1 )) . I  PHIL  IN (1) « SAVE 2  III  I  «  1351 

1  ISAVEl (3001 ),MCBLIN< 1)1 , 1  SAVE 2! 3C01I ,PHIPLN! 1)  I,  1352 

2  ( SAVE 1 16001 ) • MPCBLNI lilt! SAVE2 (6001 1 1 BIGSEM! 1 ) ) •  1353 

3  (ALIM!  lltXSMALLM  AAAAAA! 1 ) ,AC0EFF!2) ) , !BB8BBB!  1 )  ,  1354 

4  BCOEFF ( 2 ) )  1355 

TENS!X)=X*10.0  1356 

DATA!  DEC!  I  1 1 1*1, 9  )/1.0,  2.C,  3. 0,4.0,  5. 0,6,0,  7.0, 0.  0,9.0/  135  7 

DATA  LET1/4H  CB/,LET2/4H  CEG/  1358 

Kl=0  1359 

K0UNT*l  1360 

00  61  1  =  1,4  1361 

61  IEXP(I)=0  1362 

DO  60  1-1,3000  1363 

62  WLOG(I)=0.0  1364 

PHILIN!I)*0.0  1365 

BIGSEM! I  1 =0.0  1366 

EMSUPS! I  1=0.0  1367 

MDBLIN! 11=0.0  1368 

PHI PLN! I ) =0.0  1369 

60  MPDBLN! I  1=0.0  1370 

WRITE (6, 51)  1371 

READ! 5, 19) XSMALL, DEGREE,  DBL  IM  1372 

READ(5«12)FRQEXP  1373 

IF!FRQEXP.NE.0)READ(5,12)(  IEXPCK ),K*1,4)  1374 

12  F0RMATI4I 10)  1375 

19  FORMAT (3E10.0,50X)  1376 

DO  90  1=2*7  81  1377 


90  ALIMCI  )»TENSf ALIMfI-1)  )  1378 

00  70  J*l *6  1379 

00  71  1*1,9  1380 

K«KUI  1381 

71  OMEGAS! K) *DEC I I ) *ALIMf  J)  1382 

70  Kl*Kl*9  1383 

0MEGAS(55)=ALIM(  7)  .  1384 

K*l  1385 

CALL  FIGURR(55»0MEGAS,K0UNT  I  1386 

41  V*ABS(OMEGAY(K) )  1387 

X*ABS(SIGMAX(K) )  1388 

IFIY.EQ.O.O)  GO  TO  13  1389 

Y0MEGA=Y  1390 

GO  TO  20  1391 

13  IFfX.EQ.O.O)  GO  TO  10  1392 

YOPEGA*X  1393 

20  CALL  SWEEPR! YOMFGA )  1394 

CALL  FIGURR! l 7, OMEGA, MOUNT )  1395 

I F ( KGUNT .GE. 2984  ) GO  TO  40  1396 

GO  TO  BO  1397 

10  KK*K*i  1398 

Kl=KMl  1399 

IF(KI.GT.2000)  GO  TO  40  1400 

00  ll  L=KK,K|  1401 

Y= ABS (OMEGAY ( L ) )  1402 

X  =  ABS ( S IGMAXt  L ) )  1403 

IFIX.NE.C.O.OR.V.NE.O.O)  GC  TO  80  1404 

11  CONTINUE  1405 

GO  TC  40  1406 

80  IFU.GE.  ICOUNT  )  GO  TO  40  1407 

K*IU1  1408 

GO  TO  41  1409 

40  00  68  1  =  1,  3000  1410 

IFCWLOGI  n.EQ. 0.0.  AND.  PHIL  INI  I  ).EQ.O.O.ANO.PHIPLNCI).EQ.C.o.AND.  1411 

1MDBLINI I  UEC.G.O.ANO.MPDBLNI  D.EQ.C.O.AND.BIGSEM!  D.EQ.C.O.ANO.  1412 

2EHSUBSI n.EQ.J.O )  GO  TO  69  1413 

JJxMCOII  ,57)  1414 

I F ( J J.EO.O )  WR l T  E ( 6, 5 1  )  1415 

68  WRITE (6,50)  WLOG  ( I ) , R IGSEM (I ) ,MDBL INI  I) ,PHI L INfl) • EMSUBS (I) ,  1416 

1MPDBLNI I ) *  PHI PLN ( I )  1417 

69  WRITF(6,14)  1418 

CALL  NEWPLT ( DEGR EEtWLOG, PHILIN.LET2  )  1419 

IFIIEXPI  D.EO.O)  GO  TO  30  1420 

JJxlEXP(l)  1421 

00  34  I«1,JJ  1422 

34  CALL  EX PLOT (WLOG. 3000. PHIL (N*  3C00»  DEGREE.LET2 )  1423 

30  CALL  NEWPLT ( DBL IM, WLOG, MDBL  IN, LETl )  1424 

IF|  IEXP(2).EQ.O)  GO  TO  31  1425 

J J* I EXPI 2 )  1426 

00  35  I*l,JJ  1427 

35  CALL  EXPLOT  I WLOG . 3000. MOBL I N, 3000. OBL IP, LET l )  1428 

31  WRITEI6.15)  1429 

CALL  NEWPLTC  DEGREE* WLOG*  PH  IPLN  *LET2 1  1430 

IF!  IFXP(3).EC.O)  GO  TO  32  1431 

JJ* I EXPI 3 )  '1432 

00  36  1*1, JJ  1433 

36  CALL  FXPLOT (WLOG . 3000. PH  I PLN, 3000. DEGREE .LET 2 1  1434 

32  CALL  NEWPLT I DBL I M, WLOG. MPOBLN. LETl )  1435 

IF! I EXPI 4 ) • EQ.O )  GO  TO  33  1436 

J J* I EXPI 4 )  82  1437 


o  o 


c 

c 


DO  37  1=1, JJ 

37  CALL  EXPLOTIWLOG, 3000, MPDBLN, 3000,DBLIN,LETI) 

51  FORMAT ( IHl .5X.5H0MEGA, 13Xv3H/M/t I6X .6H/M/-06.  10X.3HPHI .2CX.4H/M/P, 
*10X, 7H/M/P-0B, 14X,  5HPHI—P ,/  ) 

50  FORMAT! 2X, 1PIE12.5.5X, IP1E12.5.6X. IP IE  12.5. 5X, IP  IE 12.5. 1 IX, 
*1P1E12.5,3X,IPIE12.5,9X, IP1E12.5) 

14  FORMAT! 1H1, ///////// ///////////, 25X«60HTHE  FOLLOWING  PLOTS  ARE  NCR 

*MAL  FREQUENCY  RESPONSE  DATA  ) 

15  FORMAT  (IHl, //////////////// ////,25X,  36HTHE  FOLLOWING  PLOTS  ARE  POP 

*0V  DATA  I 

33  RETURN 
END 


SUBROUTINE  F IGURR ( KONT , OMEGA. MOUNT  I 
DIMENSION  SAVEl! 100,  1001 , SAVE2! ICO. 100) 

DIMENSION  ACOEFFI IOOI.BCOEFFI 100 1 , AAAAAA! 991 .B8BBBBI99) 

DIMENSION  WLOGI 3000) .PHIL  IN ! 3000  I ,MOBL INI  30001 .PMIPLNI 30C0I , 
♦MPDBLN! 3000 I.BIGSEMI 3000 ) ,EMSUBSI3CC0I, OMEGA (KONT) 

REAL  MDBL I N, MPORLN 
INTEGER  PO INT ( 130. 100 1 
COMMON/ FREEK/ACOEFF, BCOEFF 

EQUIVALENCE! AAAAAA! 1 1 ,  ACOEFFI 2  ) ) , ! 6BBB8BI 1 ) .BCOEFF 12)  ) 

COMMON/ INF06/EMSUBS 
COMMON  SAVEl, SAVE2, POINT 

EQUIVALENCE  < WLOG(l) , SAVE  11 1 ) ) , I PHILINI l ) , SAVE2I 1 1 ) , 

1  ! SAVE  1 (3001 1 , MCBL IN! 1 ) ) • ! SAVE 2 1  3001) .PHIPLN! II), 

2  ( SAVEl (6001 ) .MPDBLN!  1 ) ), ( SAVE 2! 6001) .BIGSEMI 1) ) 

DO  21  1=1, KONT 

CALL  FIGURE (OMEGA! I )  ,  ACOEFFI 1 ) , AAAAAA, EMOFAS, PHI AS) 

CALL  FIGURE (OMEGA! I ) .PCOEFF! 1 ) ,B8BBBB, EMOFBS.PHIBS ) 
B1GEMM=ABS(EM0FAS/EM0FBS  ) 

BIGSEM! KOUNT )=BIGEMM 
WLOGI KOUNT )>OMEGAI I) 

PH  I END=PH I  AS- PHI  BS 

IF (PHIEND.LT, -180. )  PH  I EN0*PHI END* 360* 

IF(PHIENO.GT.  180.)  PHIENO*PHl END-360. 

PHILINI KOUNT )  =  PH  I  END 

PHI  END  =  PHIENO  ♦  0.0174533 

REEL  =  BIGEFM  ♦  COS  !  PHIENC  ) 

EIMAG=  BIGEMM  *  SIN  (  PHIENO  ) 

EMSUBP  *  SORT  (  (  REEL**  2  )  ♦  (  (  OMEGA  (  I  I  *  ElMAG  )  **  2  )  ) 

EMSUBS(KOUNT)=EMSUBP 

QUANZ=OMEGA( I )  *  EIMAG 

PHIPEE*ATAN2(QUANZ,REEL) 

PHIPEE  =  PHIPEE  *  57.2957795131 

IF! PHIPEE  •  LT  •  0.)  PHIPEE  =  PHIPEE  ♦  360. 

I F ( PH  I  PEE. LT. -180. )  PHIPEE=PHIPEE*360. 

IF! PHIPEE. GT.  180.)  PHlPEE*PHIPEE-360. 

PHI PLNIKOUNT ) =PH I  PEE 
DEEBEE=20.0*( ALOGIOI BIGEMM  )  I 
MDBL IN! KOUNT )*DEEBEE 
DEEBEP=20.0*( ALOGIO! EMSUBP)  ) 

MPDBLN! KOUNT )=DEEBEP 
KOUNT=KOUNT ♦! 

21  CONTINUE 
RETURN 
END 
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1438 

1439 

1440 

1441 

1442 

1443 

1444 

1445 

1446 

1447 

1448 

1449 

1450 

1451 

1452 

1453 

1454 

1455 

1456 

1457 

1458 

1459 

1460 

1461 

1462 

1463 

1464 

1465 

1466 
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1468 

1469 

1470 

1471 

1472 

1473 

1474 

1475 

1476 

1477 

1478 

1479 

1480 

1481 

1482 

1483 

1484 

1485 

1486 

1487 

1488 

1489 

1490 

1491 

1492 

1493 

1494 

1495 

1496 

1497 


n  o 


SUBROUT  INF  F  IGURE  C  X,  Y ,Z,  EMOFAS,PHI AS  I  1498 

DIMENSION  Z 199)  1499 

REELAS* V  1900 

1*2  1901 

A—1.0  1902 

21  REELAS* I(Z(I)*IX**I))*A) ♦ RE EL AS  1903 

25  l*U2  1904 

IFU.GT.99  |  CO  TO  22  1505 

A*-A  1506 

IF (Z( I )• FQ.O. ICO  TO  25  1507 

GO  TC  21  1508 

22  1*1  1509 

A*1.0  1510 

AIMAGA=0«0  1511 

24  AIMAGA*(Z(I)*f X**I)*A)*AIMAGA  1512 

26  I  *  I ♦?  1513 

IFU.GT.99  I  GO  TO  23  1514 

A=-A  1515 

IFIZm.EQ.O.  IGO  TO  26  1516 

GO  TO  24  1517 

23  EMOFAS=SCRT I ( REEL  AS* *2 1 ♦( A I MAG A* *2 1  I  1518 

PHIAS*ATAN2( Al  MAGA, REELAS )  1519 

PHI AS  *  PH I  AS *57. 2957795131  1520 

IFIPHIAS.LT. 0. )  PHIAS*PHIAS«360.  1521 

RETURN  1522 

END  1523 

C  1524 

C  1525 

SUBROUTINE  LOGGER ( XX , L IMI TW )  1526 

OIHENSION  ALIM( 7  )  »DIHS( 20)  1527 

EQUIVALENCE  I WASTE,  AXX )  1528 

COMMON/ INF02/ALIM  1529 

TENS  I X )  =  X*10. 0  1530 

OATA  IDIMSf I), 1*1,201  / 1.0, 1. 15,  1. 3, l. 5* 1. 7, 1. 95. 2.25.2. 5,2.  85 ,  1531 

*3. 15, 3.5, 3. 95, 4. 5,4. 85, 5.4, 6.05,  6. 8, 7. 7, 8*8, 9. 999999/  1532 

LIM!TW*0  1533 

L!MITX*0  1534 

LIMIT2*0  1535 

00  20  1=1,6  1536 

AA=2*f 1-1 1  1537 

20  IFIXX.LT.AL1MI !♦ I  I  .ANO.XX.GE. AL IM ( I ) )  GO  TO  21  1538 

RETURN  1539 

21  LIMI TX*TFNS( AA )  1540 

WASTE=ABS ( XX )  1541 

IF(UASTE.EQ.O.O)  RETURN  1542 

30  IFIAXX.GE.1.0.AND.AXX.LT.10.0)  GO  TC  22  1543 

IFIAXX.LT.IO.OIGO  TO  23  1544 

AXX=AXX/ 10.0  1545 

GO  TO  30  1546 

23  AXX*TENS ( AXX )  1547 

GO  TO  30  1548 

22  DO  24  1*1,19  1549 

24  IFCAXX.GE.OIMSI  II  •  AND.  AXX.LT.D  IMS!  Mill  GO  TO  25  1550 

25  LIMIT2* I  1551 

LIMITW*LIMITX4LIMIT2  1552 

RETURN  1553 

END  1554 


1555 

1556 

1557 


SUBROUTINE  BLANKR  (PO INT ,N1 ,N2 ,L IMI 1 
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o  u 


INTEGER  POINT f Nl ,N2) • BLANK,  £*iSH,0  1W8 

COMMON  / INF04/  Q, BLANK. OASH  1*59 

CAT A  BLANK/ IH  /,0ASH/1H-/,Q/IH|/  I960 

DO  91  M*1,N1  1961 

DO  91  N*1«N2  1962 

91  POINT! H,NI«BLANK  196) 

00  92  M*l,Nl  1964 

POINT IN* l I *DASH  1969 

92  POINT (M#N2 I* DASH  1966 

00  93  M*l,N2  1967 

POINT f 1,M)*Q  1968 

93  P0INT!LIM1,M)*C  1969 

RETURN  1970 

END  1971 

1972 
197  J 

SUBROUTINE  NEMPLTI X, V, Z*LET  I  1974 

INTEGER  POINT! 13C, 100) , POINTF!  120,58 1 , STAR, Q  1979 

DIMENSION  ALIM! 7 1  1976 

DIMENSION  V! 3000 )  ,Zf 3000)  1977 

DIMENSION  SAVE if 100, 1001 »SAVE2( 100#  1001  1978 

COMMON  / INF04/  0, BLANK, OASH  1979 

COMMON/ INF02/ ALIM  1980 

COMMON  SAVEl,SAVE2,P0INT  198* 

EQUIVALENCE! POINT!  1, 1), POINTF!  1,  111  1982 

DATA  ST AR/ IH*/ ,0/  1HI /  1983 

1  FORMAT! 1H11  1984 

CALL  P.LANKR!  POINTF,  120,  98,120)  1989 

DO  10  1*1,3000  1986 

IFIY!  II.EQ.O.O)  GO  TO  11  1987 

GO  TC  12  1988 

11  IFIZID.EQ.O.O)  GO  TO  10  1989 

12  CALL  LOGGER I Y f 1 1 , L IM ITN )  1990 

IFfLIMITU.EG.O)  GO  TO  10  1991 

CALL  LINEAR(X,LIM|TY,2f It)  1992 

IF! LIMITV.EC.O)  GO  TO  10  199) 

P0INTF!LIMITW,LIMITY1*STAR  1994 

10  CONTINUE  1999 

DO  40  1*20,100,20  1596 

DO  40  J*l,58  1997 

40  POINTF! I , J)*0  1598 

MRI TEf 6, 1 )  1999 

30  FORMAT! 1X,6| 1P1E 10*1, IX, 3HR/S, 6X ) , 1P1E 10, 1 )  1600 

J=X  1601 

URI TEf 6,91 )J, LET, f POINTF! I, 1 1 ,  1*1, 1201  1602 

91  FORMAT! IX, 14,  A4,1X,120A1)  1603 

WRITE f 6,201!  f  POINTF!  I,J),  I*  1,  120), J«2,57l  1604 

J*-X  1609 

WRITE!6,91)J,LET,f POINTF!  I, 98) , I *1, 1201  1606 

20  FORMAT! 10X,12CAl I  1607 

WRITE!6,30)ALIM  1608 

RETURN  1609 

END  1610 

C  1611 

C  1612 

SUBROUTINE  LINEAR! X,LIMITV, V)  161) 

DIMENSION  ENCRMT f  99)  1614 

LIMITY*0  .  1619 

FLIMIT*X*2,  1616 

DELTA*FLIMlT/58.  85  1617 


on  on 


ENCRMTI l  )*X  1618 

00  10  1*2,58  1619 

10  ENCRMTI I )*ENCRMT(  I-1I-0ELTA  1620 

ENCRMT(59)*-X  1621 

00  11  1*1,58  1622 

11  IFIY.LE.ENCRMTin.ANO.Y.GE.ENCRMTI  I  +  IHGO  TO  20  1623 

GO  TO  40  1624 

20  LIMITV*I  1625 

40  RETURN  1626 

ENO  1627 

1628 

1629 

SUBROUTINE  SWEEPR ( XOMEGA )  1630 

01MENSI0N  OMEGA! 17),SWEPERI  17)  1631 

COMMON/ INF03/0MEGA  1632 

OATA  (SHEPERl I )» 1*1, 17)  /C. Cl, 0.05,0. 1,0. 2,0. 3, 0.4, 0.5, 0.6,0. 7,  1633 

*0. 8,0.9,1.0,2.0,4.0,6.0,8.0,10.0/  1634 

00  70  K* 1, 1 7  1635 

70  OMEGA! Kl*  SWEPER I K )*XOMEGA  1636 

RETURN  1637 

END  1638 

1639 


1640 
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APPENDIX  B 


FLOW  CHART  OF  SOURCE  DECK 
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A  computer  program  to  perform  dynamic  systems  analyses  and  plot  the  results 
is  presented.  Linear  systems  and  a  large  classification  of  nonlinear  systems 
representing  engineering,  scientific,  and  economic  disciplines  can  be  modeled  to 
permit  application  of  the  computer  program.  Two  examples  are  given  to  demonstrate 
the  capabilities  of  the  analysis  tool.  The  mathematical  model  of  a  missile  guidance 
and  control  system  is  analyzed  and  a  ratio  of  polynomials  representing  the  closed 
loop  transfer  function  of  a  high  performance  model  follower  aircraft  is  evaluated. 
Linear  differential  equations  to  the  100th  order  having  real  or  complex  roots  can  be 
studied.  System  characteristic  equation  root  loci  and  system  transfer  functions 
are  plotted. 
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